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ABSTRACT 


Autonomous vehicle teams have great potential in a wide range of maritime sensing 
applications, including mine countermeasures (MCM). A key enabler for successfully 
employing autonomous vehicles in MCM missions is motion planning, a collection of 
algorithms for designing trajectories that vehicles must follow. For maximum utility, 
these algorithms must consider the capabilities and limitations of each team member. 
At a minimum, they should incorporate dynamic and operational constraints to ensure 
trajectories are feasible. Another goal is maximizing sensor performance in the presence 
of uncertainty. Optimal control provides a useful framework for solving these types of 
motion planning problems with dynamic constraints and different performance objectives, 
but they usually require numerical solutions. Recent advances in numerical methods 
have produced a general mathematical and computational framework for numerically 
solving optimal control problems with parameter uncertainty—generalized optimal control 
(GenOC)—thus making it possible to numerically solve optimal search problems with 
multiple searcher, sensor, and target models. 

In this dissertation, we use the GenOC framework to solve motion planning problems 
for different MCM search missions conducted by autonomous surface and underwater 
vehicles. Physics-based sonar detection models are developed for operationally relevant 
MCM sensors, and the resulting optimal search trajectories improve mine detection 
performance over conventional lawnmower survey patterns—especially under time or 
resource constraints. Simulation results highlight the flexibility of this approach for 
optimal motion planning and pre-mission analysis. Finally, a novel application of this 
framework is presented to address inverse problems relating search performance to sensor 
design, team composition, and mission planning for MCM CONOPS development. 
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CHAPTER 1: 
Introduction 


Over the last two decades, unmanned vehicle systems have grown steadily more capable, 
reliable, and ubiquitous. Most systems, however, are still designed to conduct specific mis¬ 
sion sets in a particular domain, with capabilities largely dependent on sensor payloads. 
As system designers increasingly turn to commercial technologies and open architectures, 
it is easier than ever for robotic systems to inter-operate. As a result, multiple dissimilar 
vehicles can be combined into a collaborative team to overcome individual vehicle lim¬ 
itations and deliver advanced capabilities—even across operating domains. Autonomous 
vehicle teams have great potential in a wide range of scientific, commercial and defense 
applications, and they are especially well-suited for remote sensing in maritime domains. 

To maximize the utility of a heterogeneous vehicle team for a given sensing mission, motion 
planning algorithms must consider the capabilities and limitations of each team member. 
At a minimum, they should incorporate dynamic and operational constraints to produce 
feasible trajectories. Optimization techniques can be used to allocate effort according to 
individual vehicles’ sensor performance. Such techniques can produce motion plans which 
are superior to conventional “lawnmower” survey patterns, which may be sub-optimal for 
certain sensors and infeasible for under-actuated vehicles to follow exactly. 

Autonomous systems must also operate with imperfect information about their environ¬ 
ment. This is particularly true in the maritime domain, where sensor accuracy usually 
depends on acoustic conditions and vehicle motion is subject to unknown disturbances at 
the water surface. In underwater search applications, the ability to detect and localize a 
target with sonar is impacted by several factors including acoustic noise, ambiguous ge¬ 
ometry, and aspect-dependence. Consequently, the performance of an autonomous system 
may depend greatly on its ability to cope with uncertainty. Motion planning algorithms 
which consider uncertainty, therefore, can increase a system’s overall robustness. 

Optimal control provides a useful mathematical framework for solving motion planning 
problems with nonlinear dynamic constraints and different performance criteria [1], [2], 
but these problems are extremely difficult to solve analytically. Pseudospectral computa- 


1 




tional methods are a powerful tool for solving these problems numerically [3], [4]. Re¬ 
cent developments in numerical methods have made it possible to explicitly incorporate 
parameter uncertainty into the objective function of an optimal control problem [5]-[10]. 
Moreover, these generalized optimal control (GenOC) problems can incorporate sensor per¬ 
formance models to produce optimal vehicle trajectories for a given sensor configuration. 
Researchers have successfully applied these methods to solve motion planning problems 
with complex, multi-agent interactions in a variety of scenarios including optimal search, 
path coverage, and force protection [11]. 

This dissertation employs the GenOC framework as a mission planning tool for different 
autonomous vehicles conducting sensing missions for mine countermeasures (MCM). In 
this chapter, we first describe the motivation for developing motion planning algorithms 
for heterogeneous autonomous vehicle teams. Next, we provide context for our technical 
approach by reviewing relevant literature in the areas of MCM, search theory and sensor 
modeling, coverage planning, and optimal control. Finally, we highlight specific contribu¬ 
tions developed during our investigation. 

1.1 Motivation 

There are a number of complex sensing missions which could utilize autonomous vehicle 
teams to deliver a mix of different capabilities. RIVERWATCH, for example, is a riverine 
environmental monitoring system comprised of an unmanned surface vessel (USV) which 
can launch and recover an unmanned aerial vehicle (UAV) to provide far-held sensing 
and improve the USV’s perceptual map of the environment for safe navigation [12]. Au¬ 
tonomous vehicle teams will play an even bigger role in missions and environments where 
safety is paramount. The U.S. Navy has embraced this vision, and has invested heavily in 
vehicle and sensor technologies for mine countermeasures (MCM). In general, MCM oper¬ 
ations are conducted in a sequence of phases, each performed by various types of vehicles 
and sensors [13]. Presently, these assets require dedicated support from manned platforms, 
but a current thrust of naval research is aimed at enabling autonomous systems to support 
other unmanned vehicles during MCM operations [14]. Since 2012, the Office of Naval 
Research (ONR) has developed a concept of operations (CONOPS) for autonomous mine 
hunting called Single Sortie Detect to Engage (SS-DTE). Eigure 1.1 provides an artist’s 
rendering of this future system in action. Under the SS-DTE concept, a USV host will 
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Figure 1.1. Concept for Single Sortie Detect to Engage for MCM operations. 
Source: [16]. 

deploy, recover, and sustain (i.e., recharge batteries, extract and process collected data, up¬ 
date mission objectives, etc.) multiple classes of autonomous underwater vehicles (AUVs) 
conducting different phases of a mine hunting operation [15]. 

Other concepts for using collaborative autonomous vehicles in MCM operations have been 
proposed. Researchers at the North Atlantic Treaty Organization (NATO) Undersea Re¬ 
search Centre (NURC), for example, have developed autonomous maneuvers for keeping 
a human-selected target in the sonar field of view of a US V to demonstrate a nascent target 
reacquisition behavior [17]. NURC researchers subsequently tested a concept whereby a 
highly capable USV detects and localizes a target in sonar imagery in order to guide an 
expendable neutralizer into the target [18]. The Naval Postgraduate School (NPS) Center 
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Figure 1.2. NPS CAVR Autonomous Surface and Underwater Vehicles. 

for Autonomous Vehicle Research (CAVR) regularly deploys fleet-representative vehicles 
and sensors in support of naval research objectives. Figure 1.2 shows two such vehicles, 
a SeaFox USV and a Remote Environmental Monitoring UnitS (REMUS) 100 AUV in 
Monterey Harbor. Motivated by their untapped potential for collaboration in MCM sensing 
missions, this dissertation develops sensor-based motion planning strategies for a hetero¬ 
geneous team comprised of these types of vehicles. 

1.2 Background 

Sensor-based motion planning means different things to different people. Choset and Bur¬ 
dick define this term as an algorithm that “incorporates sensor information, reflecting the 
current state of the environment, into a robot’s planning process” [19]. The authors dis¬ 
tinguish this approach from classical planning, which has the luxury of full knowledge 
about the world before planning begins. In this context, “sensor-based” can refer to any 
online planning algorithm. Another definition for sensor planning, or sensor path plan¬ 
ning used in the literature “refers to the problem of determining a strategy for gathering 
sensor measurements to support a sensing objective” [20]-[22]. In this dissertation, we 
extend this definition to develop motion plans which explicitly optimize the performance 
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of a given sensor on a given platform for a given mission. This work is motivated by naval 
underwater MCM missions, which are fundamentally sensing problems conducted to de¬ 
tect and classify mine targets. In this context, sensor performance refers to the probability 
of detecting a target in the available mission time. An optimal motion plan maximizes this 
performance objective. We achieve this in our formulation by minimizing its complement, 
the probability of failing to detect a target. 

1.2.1 Mine Countermeasures 

Washburn and Kress describe three methods for clearing a minefield: destruction, hunting, 
and sweeping. Destruction refers to the application of “sufficient lethal force” to destroy 
all the mines in an area. The monetary cost and severe impact to the environment of this 
approach usually make destruction a method of last resort. Therefore, MCM operations 
traditionally comprise both mine hunting and minesweeping tasks. The authors provide the 
following definitions: 

One sweeps for mines by attempting to cause the mine’s sensors to detonate the 
mine in circumstances where the detonation is harmless. If the mine is located 
by some means not involving its own sensors, and then either destroyed or 
avoided, one is instead hunting ... Sweeping and hunting are both essentially 
search problems. [23] 

The U.S. Navy (USN)’s MCM capabilities and practices are summarized in 21st Century 
U.S. Navy Mine Warfare, published by the Program Executive Office for Littoral and Mine 
Warfare. The authors state that since “minesweeping is more risky to the sweeping platform 
than mine hunting and, when completed, generally leaves behind a higher residual risk to 
vessels that transit the swept area” most MCM operational plans minimize risk by including 
both mine hunting and mine sweeping [24]. This dissertation considers only mine hunting 
operations, which are further described as follows: 

Mine hunting provides a relatively high degree of certainty that an area of con¬ 
cern is mine-free or the risk of a mine strike has been minimized. It comprises 
five steps: detection, classification, localization, identification, and neutraliza¬ 
tion. Sonars are the primary means to detect and classify mine-like contacts. 
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Identifying each contact as a mine or a “NOMBO” (Non-Mine/Mine-Like Bot¬ 
tom Object) can also be carried out by EOD divers and the Navy’s marine 
mammal systems, video cameras on mine neutralization vehicles, and laser 
systems. In this regard, advanced sonars on unmanned underwater vehicles 
offer good promise to enhance mine-hunting capabilities. [24] 

Our approach concerns the latter, since the performance of mine hunting sonars are directly 
related to the capabilities of the vehicles that deploy them. Optimizing sensor performance 
through vehicle motion is the objective of this dissertation, hence the title: Optimal Sensor- 
Based Motion Planning for Autonomous Vehicle Teams. We do not address the neutraliza¬ 
tion task, but this is often considered a special case of an identification mission whereby 
an autonomous neutralizer vehicle reacquires (and destroys) a previously-identified target. 
Mine hunting operations face several challenges due to uncertainty, namely: 

A contact that is classified as mine-like must be identified as a mine or 
NOMBO and, if a mine, rendered safe before the Navy mine countermeasures 
commander ... can declare a route or area cleared. Depending on the accuracy 
of the location of the contact, the characteristics of the bottom (e.g., smooth 
or rough), sediment type, amount of clutter, and the depth of the water, among 
other factors, the process of reacquisition and identification of each mine-like 
contact can take several hours. [24] 

We note that while the SS-DTE concept “gets the man out of the minefield” by replacing 
manned AUV support vessels with a special-purpose USV, it otherwise follows the same 
operational paradigm employed today to carry out the five steps of detection, classifica¬ 
tion, localization, identification, and neutralization. This represents a bottleneck in MCM 
operations. In fact, a Naval Research Advisory Committee report from 2000 states: 

While time lines have shortened, current sensors require repeated acquisition of 
contacts to determine whether they are mine-like or NOMBOs to make positive 
identification ... The critical issue here is that this process is presently done 
serially—one goes through the detection/classification processes by having to 
reacquire the detected objects in order to classify them, and to do so again if 
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there is a deeision to neutralize. This is not only very time eonsuming, but 
it also puts the MCM forees at risk, sinee little of the required aetion ean be 
earried out elandestinely at present. [25] 

Sinee the publieation of this report, AUVs have improved the Navy’s ability to eonduet 
elandestine aetions at greatly redueed risk to MCM personnel, but this serialized proeess 
is still time-eonsuming. We suggest that outfitting a host USV with wide-area mine detee- 
tion sonar has potential to reduee eurrent operational timelines by eondueting the first four 
steps in parallel. This dissertation develops sensor-based motion planning strategies that 
investigate these new CONOPS. 

Mine hunting tasks are often eondueted in three distinet phases: 1) wide area deteetion 
and loealization of mine-like objeets (MLOs) with long-range, low-resolution sensors; 2) 
identifieation with short-range, high-resolution sensors; and 3) neutralization [26]. Phases 
one and two are inereasingly eondueted with autonomous vehieles, while phase three is 
still eondueted primarily with explosive ordnanee disposal (EOD) divers and/or remotely- 
operated vehieles. Therefore, in this dissertation, we foeus on sensor-based motion plan¬ 
ning for phases one and two. The literature provides examples of planning methods ap- 
plieable to both phases. Phase one is typieally treated as a eoverage problem [26]-[28], 
espeeially when there is no prior information about mine loeations. Large seareh areas are 
usually subdivided into smaller regions whieh ean be eovered by a typieal AUV mission 
duration [27], or subdivided by bottom type [29], so that different sensors, traek spaeing, 
ete. ean be tailored to the loeal environment. Several eoverage path planning methods may 
be used, but vehiele paths for phase one are usually based on a deterministie eoverage 
pattern like the popular lawnmower pattern. 

Phase two ean be thought of as a targeted eoverage problem, guided by prior information 
about the expeeted target loeations. When a vehiele must visit multiple MLOs during a 
reaequire-identify (RID) mission, this is akin to solving a traveling salesman problem [28] 
to visit eaeh target loeation. Typieally, a “standard” multiple aspeet eoverage pattern eom- 
prised of parallel traeks at different headings is then exeeuted above eaeh target. This 
method ensures that high-resolution imagery is eolleeted from several different aspeet an¬ 
gles in order to aid the identifieation effort. Often, the ehoiee of sensor dietates the type of 
vehiele paths eonsidered for RID missions. Sidesean sonar, for example, does not produee 
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good imagery when the seareh platform does not follow level, straight line paths. There¬ 
fore, these planning methods do not eonsider vehicle dynamics beyond speed, minimum 
turning radius and/or the time required to make a turn. Some methods account for the dis¬ 
tance required to stabilize on the next track line after making a U-turn [30]. When planning 
lawnmower-type coverage patterns for sidescan sonar (SSS), some algorithms optimize 
over the space of track line headings [22], [31]. 

Johannsson et al. propose an algorithm for conducting both mine detection and RID mis¬ 
sions from a single AUV [32]. The planning algorithm considers the benefit of altering the 
plan to inspect a detected object immediately, or postpone identification until later, i.e., if 
the current plan will already revisit the area where the detection occurred. This approach 
is an exception to the current MCM practice of conducting phase two tasks with different 
vehicles and sensors than the ones utilized for phase one. 

The literature contains a few examples which conduct MCM tasks via heterogeneous vehi¬ 
cle teams. Sariel et al. proposes a multi-vehicle architecture in which AUVs detect MLOs 
and crawler robots revisit these locations to identify the targets with cameras. This imple¬ 
mentation provides a mechanism for resource management; the available AUVs negotiate 
amongst themselves to select from all possible waypoint locations, effectively dividing the 
search area by construction of their search paths [33]. This example does not incorporate 
realistic vehicle dynamics or sensor models, however. Bays et al. describe an automated 
scheduling system for the Navy’s SS-DTE program, which schedules a sequence of opera¬ 
tions for the team to complete, including USV transit operations to rapidly carry AUVs to 
and from their search areas [34]. 

Shafer et al. describe a developmental system comprised of US Vs and AUVs that can 
carry different sonar to cooperatively conduct phases one and two. The authors employ 
a behavior-based architecture that “performs multi-objective optimization to reconcile be¬ 
havior output” from a set of desired behaviors [35]. In this approach, each vehicle stores a 
gridded representation of the environment. When a vehicle’s sensor footprint passes over a 
cell, the “clearance level” for that cell is updated and communicated to the other vehicles in 
the team. In this manner, the authors successfully demonstrate various cooperative search 
strategies. We note that the ideas expressed by these authors are fundamentally similar to 
those presented in this dissertation. The main distinction of our approach, therefore, is the 
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direct optimization of sensor performance through vehicle motion, rather than maximizing 
sensor coverage through cooperative, emergent behaviors. 

1.2.2 Search Theory 

Search theory as a discipline in operations research grew out of work by the Antisubmarine 
Warfare Operations Research Group in World War II, first published by Koopman in 1946 
[36] and reprinted in expanded form in 1980 [37]. These principles are still widely used 
today, and are the foundation for many practical search planning methods defined in the 
International Aeronautical and Maritime Search and Rescue (lAMSAR) manual, as well 
as computer-based methods like Search and Rescue Optimal Planning System (SAROPS) 
[38]. Many of these ideas have been summarized in a report developed for the search and 
rescue (SAR) community by [39]. Stone extended these concepts to address the optimal 
allocation of search effort in [40]. In [41], Stone lists three basic elements of an optimal 
search problem: 

1. A probability distribution to quantify information about the target’s location. 

2. A detection function that relates search effort to the probability of detection. 

3. Constraints on the search effort. 

Chapter 3 of this dissertation formulates different MCM missions using various combina¬ 
tions of these elements. Stone illustrates the application of optimal search techniques for 
finding a drifting vessel lost at sea in [42], and this example problem has many of the hall¬ 
marks of MCM mission planning. Specifically, it involves sub-dividing a large search area 
into discrete cells, each with an assumed prior target distribution model, and specifying par¬ 
allel search paths according to the sensor sweep width. Detection probability as a function 
of time (for random and ideal search strategies) and as a function of parallel path spacing 
are both presented. Although searching for stationary targets (e.g., mines) is simpler than 
searching for moving targets, MCM missions are complicated by sea floor characteristics 
and clutter which can appear as mines in sonar imagery. 

Stone and Stanshine consider optimal search problems in the presence of these false tar¬ 
gets in [43], with additional theorems and solution strategies presented in [40]. In short, 
this approach assumes that search is conducted in two phases: I) broad search with a 
broad-search sensor to detect contacts, and 2) contact investigation by a sensor with a 
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much smaller sweep width to identify which contacts are actual targets. The question, 
then, is how to optimally allocate search effort between these two phases, depending on 
the contact investigation policy. Results are provided for conclusive contact investigation 
policies, which, once begun, must continue until a contact is identified as a true or false 
target. While searchers are not required to initiate a contact investigation upon each de¬ 
tection event, search strategies which implement so-called immediate contact investigation 
policies outperform those that delay this operation [40]. A recent example of this approach, 
applied to a military surveillance operation, utilizes a verification team to investigate con¬ 
tacts reported by an imperfect sensor in order to find a target in minimum time [44]. 

Stone’s solution strategies for problems with false targets mirror MCM search and RID 
missions. We remark, however, that current operational practice does not conduct these 
missions concurrently. Nevertheless, Stone’s results suggest that efficiencies can be real¬ 
ized by cooperative, heterogeneous vehicle teams that can implement immediate contact 
investigation. Since autonomous systems are reaching a level of maturity that makes such 
cooperation possible, this is an area worthy of future development. While this disserta¬ 
tion describes a method for generating optimal sensor-based search trajectories for a team 
of dissimilar vehicles and sensors, the method does not yet accommodate motion plans 
to conduct MCM search and RID missions concurrently. This capability would allow mis¬ 
sion planners to address false targets through immediate contact investigation in the manner 
of [40]. At present, however, we operate under the prevailing operational assumption that 
separate RID missions are required to positively identify mine targets from a list of MLOs 
detected by an initial survey. We then utilize the GenOC solution framework to generate 
motion plans for both types of MCM missions. 

The field of search theory is vast and continuously evolving. Benkoski provides a “com¬ 
prehensive review of the existing published search theory literature,” with specific applica¬ 
tions divided according to “one-sided search” for non-evading targets, and “search games” 
whereby targets may work to thwart the search operation [45]. More recently, Chung et 
al. surveyed the results in this field with particular relevance to search operations by mo¬ 
bile robots. These authors generate a “partial taxonomy of the parameter space for search 
and pursuit-evasion models” [46]. This taxonomy includes branches for searcher, target, 
and environment models. The searcher branch incorporates models that quantify searcher 
motion, sensor detection, and the number and type of search agents available. The target 
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branch categorizes applications by the number and type of targets being hunted, whether 
targets are stationary or mobile, and if mobile, whether their motion is adversarial. Envi¬ 
ronmental models, meanwhile, are subdivided into continuous and discrete representations. 
Under this taxonomy, the search problems considered in this dissertation include: 

• one or more heterogeneous vehicles and sensors, conducting a one-sided search for 

• one or more stationary targets, with known or unknown (uniform) prior distributions, 

• in continuous time and space. 

Historically, these problems derive from Koopman’s probabilistic search results and de¬ 
fine measures of performance based on expected value [46]. This dissertation uses the 
same theoretical approach, developing an objective function in Chapter 3 to quantify the 
expected probability of failing to detect a target. Optimal search trajectories that minimize 
this objective function in continuous time and space can be generated via numerical ap¬ 
proximation techniques in the GenOC framework. In fact. Stone et al. refer to “optimal 
search in continuous time and space subject to constraints” as uncertain optimal control 
problems [47]. These techniques distinguish our approach from optimal search strategies 
which rely on discrete representations of the environment, searchers, and/or sensors. 

Hollinger et al. propose methods for distributed, multi-vehicle search planning that al¬ 
low searchers to share and fuse data about a target’s position and improve team perfor¬ 
mance [48], [49]. Realistic communications constraints are enforced by a physics-based 
model of the acoustic communications channel, but no motion constraints are imposed on 
searcher motion, detection is simulated by ideal sensor models, and search plans consist 
of transitions between discrete grid cells. Another information-based, discrete search for¬ 
mulation is presented by Bay log and Wettergren in [50]. The authors model the intake 
of sensor information during a search operation as an information flow through a com¬ 
munications channel. This approach replaces traditional sensor models with information 
measures that determine how much information can be gained from visiting discrete cells 
in the search space. The resulting optimal search plan comprises an optimal cell visitation 
sequence. While this framework can incorporate true target detections and false positives 
as a function of the channel characteristics in each cell, this approach does not address 
searcher dynamics. 
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1.2.3 Sensor Performance Modeling 

Our definition for sensor-based motion planning requires a model of the sensor whose per¬ 
formance we seek to optimize. For underwater applications, sonar is the sensor of choice. 
The MCM missions discussed thus far employ several different types of sonar, and mea¬ 
sures of performance might therefore depend upon a specific mission objective. Operators 
planning a wide area survey with a long-range sonar may be concerned only with detect¬ 
ing MLOs. Meanwhile, operators planning a RID mission with high-resolution SSS or 
multibeam sonar may be concerned only with achieving complete coverage of a designated 
area of the sea floor. Even though several algorithms have been proposed that guarantee 
complete coverage, this objective might not ensure even coverage, e.g., due to localized 
environmental conditions. Therefore, we seek an objective function that can be applied to 
different sensing objectives, yet still model realistic performance. In this dissertation, we 
address phase one and two MCM missions as optimal search problems, and these prob¬ 
lems typically require some type of detection model. Therefore, we make the following 
assumptions with regard to mission performance: 

1. Phase one succeeds if we detect a target with a long-range sensor. 

2. Phase two succeeds if we detect a previously detected target with a high-resolution 
sensor. 

Many of the detection models developed for search theory by Koopman [36] are still in 
wide use today. Washburn provides an excellent overview in [51], including models based 
on definite range, lateral range curves, and the inverse cube law, as well as detection rates 
based on “glimpses” with a physical sonar or radar model. The models derived in Chap¬ 
ter 2 of this disseration are based on signal excess and utilize the latter approach. Hansen 
provides a useful introduction to the principles of underwater sound for different sonar 
applications in [52]. This information is helpful in understanding the concept of signal 
excess. For the interested reader, a number of detailed texts on this subject are available, 
including [53]-[56], to name a few. 

The sonar detection models developed in this dissertation were designed to rapidly calcu¬ 
late signal excess within a trajectory optimization routine, so complex environmental and 
reverberation effects are neglected. Other sonar performance models, including the U.S. 
Navy’s Comprehensive Acoustic System Simulation (CASS), can model these effects via 
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different mathematieal models tailored to a given operational environment. It has been used 
to analyze the effeets of uncertain environmental data on mine detection performance [57], 
[58], but this model is too computationally-intensive for our application. 

There are a number of existing sonar perfomance models to choose from (Etter lists 
twenty six! [56]), but many do not lend themselves to motion planning, particularly if 
environmental effects are desired. Once exception is Planning Aid for Tasking Hetero¬ 
geneous Assets (PATHA), which uses Extensible Performance and Evaluation Suite for 
Sonar (ESPRESSO), NATO’s standard minehunting sonar performance prediction tool. 
This planning tool optimally assigns vehicles to straight line track segments according 
to their sensing capabilities and the local sea floor characteristics [59]. PATHA estimates 
detection probabilities by generating lateral range curves for a given sonar, target, seabed 
type, and vehicle platform; it does not generate optimal trajectories for these vehicles to 
follow, but simply specifies a location for their desired track lines. 

Einally, some applications analyze sonar performance based on the imagery they produce, 
and have developed models that simulate image construction by SSS and high-resolution 
forward-looking sonar (EES) [60], including modern imaging sonars like the BlueView 
P900 EES [61], [62]. These methods can be useful for simulating the performance of 
computer-aided detection and classification algorithms on synthetic sonar imagery, accel¬ 
erating the development of new capabilities when actual sonar data is scarce. 

1.2.4 Coverage 

Eor many robotic sensing applications, the stated objective is complete coverage of an 
operating area with a sensor. A typical example is conducting a survey with an AUV to map 
a region of the sea floor with high-resolution sonar imagery. In this dissertation, however, 
the objective of sensor-based motion planning is to optimize the detection performance of 
a given sensor. The difference is nuanced, but implies that a sensor’s performance can vary 
along a path that was designed to obtain complete coverage. This is realistic, as we know 
that factors like vehicle motion or localized environmental conditions at points along the 
robot’s path can adversely impact its sensor performance. 

Cai and Eerrari propose a sensor planning algorithm for classifying multiple fixed targets 
by strict placement of a robotic sensor. The algorithm explicitly accounts for sensor place- 
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ment by defining a subset of the robot’s eonfiguration spaee from whieh the sensor’s field of 
view (FOV), rigidly attached to the robot, will intersect each target. It then uses an approx¬ 
imate cell-decomposition method to find a path that maximizes an information-theoretic 
metric. The proposed application is similar to an MCM RID mission: a ground robot is 
tasked with revisiting targets previously detected by an aircraft’s infrared sensor, and de¬ 
ploying a ground penetrating radar to positively identify each target [20]. A weakness of 
this algorithm is that near-perfect information about the robot and its environment are re¬ 
quired. While this situation might be possible for ground robots, it is certainly not the case 
for underwater vehicles. Nevertheless, this algorithm is representative of motion planning 
based on a sensor’s FOV, or footprint; but not necessarily on its performance. The sonar 
detection models developed in this dissertation incorporate three-dimensional FOVs based 
on realistic sonar design criteria, however these FOV limits are reflected primarily through 
the sonar’s detection performance. 

Coverage planning is another method which has historically focused more on sensor place¬ 
ment than sensor performance. A number of coverage planning algorithms for robotics 
exist. In his survey of this literature, Choset describes “applications such as floor clean¬ 
ing, lawn mowing, mine hunting, harvesting, painting, etc. [that] require a coverage path 
planning algorithm, which as its name suggests, specifically emphasizes the space swept 
out by the robot’s sensor” [63]. Historically, coverage planners first decompose a two- 
dimensional planar region into cells which are easy to cover individually, then achieve 
provably complete coverage by planning a path that visits each region. A famous ex¬ 
ample of this approach is the boustrophodon exact cellular decomposition by Choset and 
Pignon [64], which allows each cell to be covered by simple back-and-forth motions in an 
“ox-plowing” or lawnmower pattern. Latimer et al. have extended this approach to multiple 
vehicles [65]. Similar approaches have produced coverage path planners which optimize 
the robot’s path length or minimize the number of turns, e.g., [66]. A recent survey of 
this literature in [67] reports on efforts to develop path planning approaches to cover three- 
dimensional surfaces such as urban structures and underwater bathymetry [68]. 
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Coverage planning algorithms are used extensively to plan underwater surveys with imag¬ 
ing sonar. As Williams explains: 

The path-planning approach currently used in practice for underwater mine 
detection operations is to design a simple ladder survey with equidistant tracks 
over the mission area. The adherence to traversing parallel tracks in MCM 
operations is partly because the collected raw data is subsequently processed 
into imagery (e.g., synthetic aperture sonar (SAS) imagery), for which such 
data is preferable. Thus, the crux of the path-planning task becomes how to 
design the spacing of (parallel) tracks. [69] 

Conceding that parallel survey tracks are often required to satisfy sonar imagery require¬ 
ments, Williams stresses that historical coverage planning algorithms which produce equal 
track spacing are sub-optimal for MCM operations. Specifically, they do not account for 
sensor performance along a given track line. Mine detection performance is usually a func¬ 
tion of highlight-to-shadow patterns in sonar images that indicate an object proud of the sea 
floor. As these patterns depend on range, geometry, and seabed type, Williams proposes 
an algorithm to determine the optimal track spacing based on these characteristics, in con¬ 
junction with SAS data. The algorithm outperforms uniform spacing, as measured by the 
proportion of search area with detection probability above a given threshold [69]. 

Pauli et al. develop sensor-driven online coverage algorithms designed to optimize the in¬ 
formation collected by a sensor. These algorithms utilize information gain from sidescan 
sonar, modeled by classification confidence [21] or detection probability [22] (both ob¬ 
tained from lateral range curves of the environment computed by ESPRESSO). While these 
algorithms incorporate sonar performance data in the manner of [69], they have the added 
advantage of adapting their patterns to align survey tracks on headings which maximize 
their objective function. Planning is accomplished on a discrete hexagonal decomposition. 
Another multi-objective coverage planning method which accounts for sensor detection 
performance via lateral range curves is described by Abreu and Matos [26]. The authors 
claim that this approach, which uses evolutionary algorithms to optimize the AUV’s sur¬ 
vey path, obtains higher detection probabilities than the optimal track spacing algorithm 
suggested by [69]. 
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Yetkin et al. propose a method for estimating the number of targets in each cell of a dis¬ 
crete grid. This algorithm assumes that probabilistic sensor performance corresponds to 
the environment in each cell, and these environments are drawn stochastically from a fi¬ 
nite set of possible environments. The authors propose a “decision-theoretic cost function 
for realistic search problems that include the effects of false alarms, missed detections, 
uncertainty in the environment, and any prior information” about the number of objects 
at a location [70]. As with the other algorithms described above, sidescan sonar motion 
constraints are imposed, so this algorithm determines the best rows to traverse in the time 
available. After sampling the first cell in its current plan, it adaptively plans the first N steps 
of a new optimal path in a receding horizon fashion. Simulations reveal that this adaptive 
approach outperforms a traditional lawnmower pattern. While no other MCM planning or 
coverage method I encountered in the literature incorporated all of these factors on search 
performance, the current implementation of this approach lacks realistic sensor models and 
vehicle dynamics. 

All of these approaches implement an exhaustive area coverage pattern by a set of straight- 
line survey tracks. In de-mining applications, Acar et al. observe that “exhaustive coverage 
is the best strategy when the robot has unlimited time (less than what is required by a ran¬ 
dom strategy) and no a priori information” [71]. When time or resource constraints prevent 
complete coverage, however, these authors recommend a probabilistic planning method 
to discover and/or exploit any prior information about the minefield itself. The rationale 
for this approach is that minefields typically have a regular pattern arising from traditional 
mine-laying methods. Once this pattern is discovered, a robot can conduct an opportunistic 
search of expected mine locations. Kim and Healey reach a similar conclusion in [72], and 
provide detailed simulation results comparing random and deterministic search strategies 
for mine clearance operations by robots with different capabilities. 

While the sensor-based motion planning algorithm in this dissertation has commonalities 
with coverage planning methods that attempt to maximize sonar performance, our approach 
does not rely on straight-line survey tracks. Indeed, the GenOC framework can accommo¬ 
date complete six-degrees of freedom (DOF) models of search vehicle dynamics to generate 
feasible vehicle trajectories that optimize detection performance. 
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1.2.5 Motion Planning via Optimal Control 

LaValle states that “a fundamental need in robotics is to have algorithms that convert high- 
level specifications of tasks from humans into low-level descriptions of how to move. The 
terms motion planning and trajectory planning are often used for these kinds of problems” 
[73]. Under this broad definition, the field of robot motion planning encompasses a number 
of widely different applications. Like search theory, the literature dedicated to this field is 
already immense. Unlike search theory, however, this field is growing at an exponentially 
faster rate. Many texts like [73] provide a good introductory overview of this subject, so 
we restrict our discussion to motion planning algorithms based on optimal control. 

Hurni presents useful definitions that help distinguish between path planning, motion plan¬ 
ning, and trajectory planning in [74]. Under these definitions, path planning finds feasible 
paths from one robot configuration to another; motion planning computes these paths, pa¬ 
rameterized by time; and trajectory planning calculates the complete evolution of a robot’s 
state and actuator control variables as a function of time. Therefore, the techniques de¬ 
scribed in this dissertation implement trajectory planning; they generate vehicle state and 
control trajectories that optimize a search objective. However, we continue to refer to the 
outputs of this process as motion plans in the context of MCM mission planning. 

Hurni also compares several path/trajectory planning methods applied to the solution of 
a benchmark obstacle avoidance problem. Multiple bug algorithms are implemented, in 
addition to various artificial potential field, roadmap, cellular decomposition, and optimal 
control-based planning methods. When evaluated against several criteria including opti¬ 
mality, feasibility, computational complexity, and portability, Hurni suggests that optimal 
control trajectory planning methods (solved via numeric optimization) outperform all other 
methods tested against his benchmark problem [74]. 

A standard optimal control problem seeks a control trajectory u{t) that forces a dy¬ 
namic system x{t) = f {x(t),u(t),t) from its initial state x(to) = xq to a desired fi¬ 
nal state x(tf) = Xf while minimizing some cost functional. A Bolza cost functional 
J[x(-),u(-),to,tf] comprises an end-point cost E (T(to)A(t/),toT/) and a running cost 
R {x{t),u{t),t) accumulated over the time domain t 6 \tQ,tf] . In its most general form, 
this problem is [75]: 
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Standard Optimal Control Problem 


minimize A 

u(t) 

subject to e(^x(to),x(tf),to,tf^ 

h (x(t),u(t)) 

For all but the simplest optimal control problems, solutions must be obtained by computa¬ 
tional methods. “Pseudospectral optimal control theory, a term coined by Ross, is a joint 
theoretical-computational framework” for solving constrained, nonlinear optimal control 
problems [4]. Ross and Karpenko review the theoretical foundations for this framework, 
discuss the ramifications for selecting a pseudospectral method, and describe practical im¬ 
plementations which utilize this approach [4]. At its heart, this solution framework trans¬ 
forms an optimal control problem from the physical domain to a computational domain 
where it can be approximated by a grid of interpolation nodes (and weighting functions) 
that have been carefully selected for convergence. Additional theoretical results are pro¬ 
vided in [3], [75]-[77]. 

Recent advances in numerical methods have made it possible to explicitly incorporate pa¬ 
rameter uncertainty into the objective function of an optimal control problem. This situa¬ 
tion arises when conducting an optimal search for stationary targets at unknown locations, 
or mobile targets whose motion can be conditionally-determined by an uncertain param¬ 
eter. Foraker develops a framework for solving such problems in continuous time and 
space as “generalized optimal control problems,” and proves that his proposed “discretiza¬ 
tion schemes are consistent approximations” to the original problem [78]. Foraker et al. 
expand upon these results in [6], [7]. The objective functional for this optimal search prob¬ 
lem “involves an expectation of a Bolza-type cost functional over a space of stochastic 
parameters” co £ Q <z described by the continuous probability density function (PDF) 
0(m) : R. This yields the generalized optimal control (GenOC) problem [5]: 


= E {x{tQ),x{tf),tQ,tf) + / R {x{t),u{t),t) 

Jta 

= x(t) (the state dynamics) 

—^ 

= 0 (the end-point conditions) (1.1) 

< 0 (the state-control path constraints) 
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Generalized Optimal Control (GenOC) Problem 


minimize 

u(t) 

J[x(-),u(-),Tf]= f E 
Jsi 

{^x(Tf),io^ -1- G |y R {x(t),u(t),t,co) Jrj 

1 (f)(io)dd) 


f{x(t),u(t)) = x(t) 

(the state dynamics) 


subject to 

to 

II 

t° 

1 

o 

tK 

(the initial condition) 

(1.2) 


h {u(t)) < 0 

(the control constraint) 



The time domain t 6 [OjT/-] is used to indicate a specified mission duration Tf, and G(-) 
is a single-valued function of the integrated running cost R Subsequent work 

by Phelps et al. have produced a computational algorithm (with consistency proofs) for 
discretizing the parameter space to produce a family of standard optimal control problems 
that can be “solved using pseudospectral discretization in the time domain” [5]. Theoretical 
foundations for this solution approach are summarized in [47]. Walton further generalizes 
this framework to problems “with parameter uncertainty in both the cost and dynamics of 
optimal control problems,” and analyzes the consistency of states and costates propagated 
by control trajectories [79]. To support these efforts, Walton developed a flexible software 
package for specifying and solving these problems via MATLAB. This software allows 
users to specify different discretization methods, e.g., Euler or pseudospectral, in the time 
and parameter domains, and “automates the creation of features which are crucial to effi¬ 
cient implementation such as gradient information and sparsity patterns” [79]. All of the 
simulation results in this dissertation were produced using Walton’s software package. 

Finally, we note that Ross et al. define a similar framework for solving optimal control 
problems with uncertain, or “tychastic” parameters. The authors refer to these problems 
as Riemann-Stieltjes optimal control problems, because their cost functional is given by 
a Riemann-Stieltjes “functional of functionals” [80]. When applied to an optimal search 
problem, this formulation produces the same objective function as Equation (1.2), which 
can be solved via the same computational methods described in [5]. A recent implementa¬ 
tion of this framework to solve a problem with uncertain dynamic constraints is described 
in [9]. 
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1.3 Contributions 

Recent theoretical results in [5], [78], [81] have produced a general mathematical and com¬ 
putational framework for solving optimal control problems with parameter uncertainty— 
generalized optimal control (GenOC). Leveraging a numerical toolbox developed in [79], 
we formulate and solve several motion planning problems for MCM operations by au¬ 
tonomous vehicles. The application of GenOC to this important naval mission set is a 
novel contribution of this dissertation. In addition, this dissertation: 

1. Develops physics-based sonar detection models to make the GenOC solution frame¬ 
work operationally relevant to the vehicles and sensors currently employed for MCM. 

2. Demonstrates, through specific examples, the utility of this method for both optimal 
trajectory planning and pre-operational mission analysis. 

3. Highlights how the proposed approach improves mine detection performance vs. 
conventional planning methods, especially under time and resource constraints. 

4. Illustrates how optimal solutions obtained from this framework can be used to solve 
inverse problems related to a specific sonar design or vehicle team configuration. 

These features make GenOC an attractive technique for investigating new MCM concepts 
of operations (CONOPS). Each contribution is presented sequentially by chapter^, as sum¬ 
marized below. 

Chapter 2 derives probabilistic detection models for various mine detection and identifica¬ 
tion sonars carried by MCM vehicles such as the Mk 18 Mod 2 Kingfish and Mk 18 Mod 1 
Swordfish AUVs. These models utilize existing techniques from sonar engineering and op¬ 
erations research for predicting sonar performance based on signal excess (computed from 
the physics-based sonar equations for a given sonar design) and detection rate. Foraker 
implemented simple models of this form to solve optimal search problems in the GenOC 
framework [78]. Walton extended this approach to applications beyond optimal search by 
replacing detection rates with attrition rate functions, deliberately shaped to model desired 
multi-agent interactions [79]. We present a method for shaping each sonar’s detection rate 
to accurately reflect its physical field of view in three dimensions. This contribution ex¬ 
plicitly relates a sonar’s performance to the trajectory executed by its search vehicle. Long 

'Portions of Chapters 2 through 4 were presented at OCEANS 2016 MTS/IEEE Monterey on September 
21, 2016 and will be published in ISBN 978-1-5090-1527-6 (IEEE Catalog Number CEP160CE-P0D). 
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acknowledged in the literature as an important aspect of MCM mission planning, this fact 
has historically been addressed only by heuristic search patterns. 

Chapter 3 presents a flexible method for formulating different MCM operations as opti¬ 
mal search problems readily solved by the GenOC framework. The method’s flexibility 
derives from its modular approach. Depending on the prior information available, differ¬ 
ent target distribution models can be constructed to produce motion plans ranging from 
wide area mine detection surveys to target identification missions. Similarly, models for 
different search vehicles and sensors can be combined to identify the most effective multi¬ 
vehicle team for a given MCM scenario. Chapter 4 demonstrates the ability of the GenOC 
framework to solve sensor-based motion planning problems for different single- and multi¬ 
vehicle MCM missions. A single-vehicle example illustrates how optimal trajectories out¬ 
perform conventional coverage patterns when the time available for conducting a search is 
limited. 

Finally, Chapter 5 proposes a novel use of the GenOC framework to investigate inverse 
problems concerning the best vehicle or sensor configuration to use in a given MCM sce¬ 
nario. The ability to rapidly solve multiple optimal search problems makes Monte Carlo 
analysis possible, providing engineering insights about the search assets themselves. This 
capability stems from the choice of sonar models derived in Chapter 2, which reflect spe¬ 
cific sonar design parameters. Hence, it is possible to analyze the impact of parame¬ 
ter changes on optimal search performance and recommend improvements that yield the 
biggest potential payoff. This dissertation provides examples of several inverse problems 
and highlights an added benefit of this approach. Namely, since each simulation produces 
a set of optimal search trajectories for the system under test, inverse problem analyses can 
be used to construct a library of optimal motion plans for a wide array of MCM missions. 
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CHAPTER 2: 

Sonar Detection Models 


In this chapter, we develop probabilistic models for two types of sensors routinely deployed 
during mine countermeasures (MCM) operations: forward-looking sonar (FLS) and sides- 
can sonar (SSS). While these sonar systems can differ widely according to their intended 
application, they share several common characteristics. First, both FLS and SSS are exam¬ 
ples of active sonar. That is, they transmit acoustic signals into the water and process the 
echoes reflected from objects in the environment to detect their presence. For a given sonar, 
this process occurs at an average rate, so detection performance depends on time. Second, 
active sonar systems employ transmit/receive arrays of transducers to improve detection 
performance in a desired direction and often add acoustic baffling to reject echoes from 
unwanted directions. A sonar design’s array geometry therefore produces an effective field 
of view (FOV) within which targets can be reliably detected, so detection performance 
also depends on a sonar’s orientation relative to targets in the environment. Last, since 
these sonar systems are rigidly-mounted onto a vehicle platform, the sonar’s orientation 
ultimately depends on that vehicle’s trajectory through the search area. This trajectory de¬ 
fines the position, orientation, and velocity of the vehicle as a function of time. Since these 
quantities are constrained by the vehicle’s equations of motion, we note that a sonar’s over¬ 
all detection performance is a function not only of its design parameters but also its vehicle 
platform dynamics. 

Assuming that detection performance defines a sonar’s effectiveness for a given mission, 
this metric can be generalized to any sonar. That is, a mission to detect and localize mines 
with a long-range, low-resolution sonar has the same objective as a mission to reacquire 
and identify these mines with a high-resolution sonar. We assume that detection with such 
a sonar is sufficient for successful identification to occur. Under these assumptions, sensor- 
based motion planning algorithms for mine countermeasures (MCM) should employ sonar 
detection models with the following characteristics: 

• Detection probability reflects an actual sonar’s dependence on array design, vehicle 
dynamics, and three-dimensional search geometry. In this way simulation can serve 
as a powerful tool to evaluate the effectiveness of different sonar designs deployed 
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from various vehicle platforms. 

• Simulated detection performance agrees with the expected/observed performance of 
actual sonar. This can be assessed by model verification and validation. 

• Detection functions permit rapid calculations within numeric optimization routines. 
This requires trade-offs between accuracy and execution speed. Smooth, differen¬ 
tiable functions with analytic gradients, for example, significantly reduce solution 
times when using gradient-based optimization. 

Many sensor models commonly used in search theory are chosen for their computational 
simplicity and do not satisfy all of these desired characteristics. Most ignore three- 
dimensional geometry, for example, but this can greatly impact detection performance 
when searching for mines on the sea floor with a surface craft’s FLS (Figure 2.1) or an 
underwater vehicle’s SSS. It is also well-known that SSS can not detect targets located 
directly beneath a vehicle’s path of travel, in the so-called near-nadir region (Figure 2.2). 
For this reason, overlapping sensor swaths are required to obtain complete coverage with 
this sensor [82]. 

Definite range models, or “cookie cutter” sensors, simply assume that detection is certain 
within a fixed range of the sensor and impossible outside it. Washburn and Kress note 
the appeal of such a model for analysis, but acknowledge that “attempts to forecast fixed 
ranges in the real world are often disappointing,” remarking that “forecast detection ranges 
for sonars are notoriously subject to error—it is not uncommon to be off by a factor of two 
or more” [23]. 

One alternative to definite range models are so-called lateral range curves. This approach 
graphs a sensor’s detection probability as a function of lateral range, defined as the distance 
from a searcher’s straight line track at its point of closest approach to a target. The area un¬ 
der a sensor’s lateral range curve defines its sweep width, a measure of sensor effectiveness 
used when planning to search an area with evenly-spaced track lines, e.g., with a lawn- 
mower coverage pattern [83]. Lateral range curves can be derived analytically, assuming 
detection rate is proportional to range via an inverse cube law, or derived empirically via 
repeated experiments [51], [84]. 

Figure 2.3 depicts the lateral range curve and corresponding sweep width for a typical 
sensor and for the special case of a definite range sensor with radius R overlaid in green. 
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Figure 2.1. Problem geometry for a surface craft with bow-mounted forward- 
looking sonar. 

We note that in general, both sensor models produee maximum deteetion probabilities at 
target ranges approaehing zero. While appropriate for optieal sensors (this was one of 
the original motivations for deriving the inverse eube law in World War II), these models 
usually require modifieations to aeeurately simulate sensors like SSS. 

One example of a sensor profile used to approximate the expeeted eoverage pattern of SSS 
(Figure 2.4) is presented in [30]. This model modifies a “eookie eutter” sensor by adding a 
“blind zone” of zero deteetion probability in the near-nadir region below the vehiele. The 
resulting sensor profile is shown in Figure 2.5. This model is well-suited for this speeifie 
applieation, namely finding an optimal traek line loeation that maximizes the probability 
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Figure 2.2. Problem geometry for an underwater vehicle with sidescan sonar, 
showing negligible sensor coverage directly below the vehicle. 



searcher track 


Figure 2.3. Lateral range curves and sweep widths for a typical sensor (black) 
and a definite range sensor (green). Adapted from [23]. 
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Figure 2.4. Sidescan sonar image 
showing near-nadir gap. Source: [85]. 



Figure 2.5. Sensor profile modified with 
near-nadir gap. Adapted from [30]. 


of detecting targets for a given track line direction, but it does not generalize well to other 
applications. Furthermore, this sensor model does not depend explicitly on vehicle dynam¬ 
ics; the knowledge that SSS performs better when its vehicle platform follows straight line 
tracks is implicit in the problem formulation, which considers only straight path segments. 


An engineering-based approach to modeling sensors like radar and sonar calculates “sig¬ 
nal excess” from physical models of the sensor and its operating environment to determine 
when detection is possible [51], [84]. Moreover, the Poisson Scan model described in [51] 
and [86] can be used to derive a sensor’s detection rate. Sensor models of this form are 
used to solve optimal search problems in [7], [11], [78], although these models implement 
a much simpler approximation of an actual signal excess equation. Nevertheless, [79] de¬ 
scribes how models based on rate functions can be calibrated to “shape” their performance 
and solve a wide variety of problems, highlighting the flexibility of this modeling approach. 
The following sections present a detailed derivation of the sonar detection models used in 
the remainder of this dissertation. The main technical contribution of this chapter is a signal 
excess model that incorporates sonar design parameters and three-dimensional geometery 
to compute detection probability as a function of a vehicle’s trajectory. 
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2.1 Signal Excess 

The signal excess model of sonar detection simulates the conditions under which an ac¬ 
tive or passive sonar system can detect an underwater object, based on well-known sonar 
equations. First proposed in [36], signal excess is still widely used in many sonar per¬ 
formance models today, including the U.S. Navy’s CASS, a Navy Standard maintained 
by its Ocean and Atmospheric Master Library (OAML) [87], [88]. CASS is a modular 
software framework which can interact with multiple different mathematical models, e.g., 
the Navy’s Gaussian Ray Bundle (GRAB) eigenray sound propagation model, to compute 
individual terms in the sonar equation. For this reason, CASS is an example of a model 
operating system (MOS), which can solve complex sonar performance modeling problems 
by decomposing them into smaller sub-problems, each tailored to a specific sonar con¬ 
figuration or operational environment. The flexibility of a MOS can be invaluable when 
designing a new sonar system or developing a tactical decision aid (TDA) to optimize a 
sonar’s performance in a given environment [56]. Their complexity and large computa¬ 
tional runtime, however, make these systems unsuitable for sensor-based motion planning 
algorithms. Anecdotally, when CASS “is used as the acoustic calculation engine ... com¬ 
putation of signal excess in support of a complex multistatic active sonar analysis task can 
take days,” with the “vast majority of runtime” devoted to computing eigenrays connecting 
a sonar transmitter, sound-scattering features, and a sonar receiver [89]. 

Our implementation, therefore, makes simplifying assumptions in order to rapidly compute 
signal excess for an active sonar attached to a moving vehicle platform: 

• We assume detection performance is limited only by acoustic background noise and 
neglect reverberation due to the sonar’s backscattered acoustic energy. Reverberation 
is a complex function of time, range, and the environment itself (e.g., the roughness 
of the sea bed) [54]. By using the noise-limited form of the active sonar equation, 
a constant parameter called the figure of merit (FOM) can be computed for a given 
sonar design, facilitating qualitative comparisons of detection performance in a given 
environment. 

• We assume that the environment is homogeneous, with a flat bottom type and con¬ 
stant water depth. We further assume a constant sound velocity profile, even though 
the speed of sound actually varies as a function of temperature, salinity, and depth. 
These assumptions allow us to forgo computationally expensive ray-tracing calcula- 
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tions. 

• We ignore multipath propagation effects, as MCM sonars typically operate at higher 
frequencies and relatively short ranges, e.g., hundreds of meters as opposed to tens 
of kilometers for a submarine sonar system. 

The signal excess model assumes that detections can only occur when the acoustic energy 
transmitted by a sonar is sufficient to overcome the two-way propagation losses in the 
environment, and the received signal reflected by a target exceeds a detection threshold 
relative to the prevailing background noise. This signal excess can be computed using well- 
known sonar equations, “simple algebraic expressions used to quantify various aspects of 
sonar performance” [84]. Because many terms in the sonar equations represent ratios of 
measured physical quantities to standard reference values, it is convenient to express them 
in units of decibels (dB), defined as lOlogjo {Ratio). For example, the reference intensity 
(power per unit area) for sonar is defined as the intensity of a planar sound wave with 
root mean square (RMS) pressure of one micropascal (/iPa). Table 11.1 in [56] defines 
several sonar parameters and their corresponding reference values. A typical signal excess 
equation for an active sonar operating against a noise background is given in [84]: 

SE = SL-IPL + TS - {N - AG) - DT (2.1) 

where SE = signal excess, 

SL = source level, 

PL = one-way propagation loss, 

TS = target strength, 

N = omni-directional sonar self-noise, 

AG = array gain, 

DT = detection threshold. 

An equivalent expression that includes signal processing terms for computing the detection 
threshold with frequency modulated (FM) or continuous wave (CW) active pulse types can 
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be found in [54]: 


SE = SL- 2PL + TS-(N-DI + lOlogio B) - DT 
where DI = directivity index, 

DT = 5 logio d-\0 logio BT -5 logjo n, 

B = pulse bandwidth [Hz], 

T = pulse duration [s], 
d = detection index, 

n = number of pings used in detection decisions. 


( 2 . 2 ) 


Recall that each term in Equation (2.2) is expressed in dB units, unless otherwise specified. 
We will use this form of the sonar equation in the sections to follow, which derive the 
physics-based probabilistic sensor models at the heart of our motion planning algorithm. 
The first step combines terms related to a specific sonar design into a figure of merit (FOM). 
This metric simplifies performance calculations during trajectory optimization, but permits 
the direct comparison of different sonar models in a given scenario. 


2.2 Figure of Merit 

For sonar performance analysis, the terms of the signal excess equation are often combined 
into the following form: 


SE{t) = EOM - PL (2.3) 

where EOM = figure of merit, 

PL = one-way propagation loss, 

D{t) = distance to the target. 

In this form, the one-way propagation loss is a function of the distance between a sonar 
and a target location m at each moment in time: D(t) = ||m - £(0||- In our MCM search 
problem, the target location m is uncertain but characterized by the probability density 
function (p{co) : i-^ R (see Section 3.3). Figure of Merit (FOM), on the other hand, is 

a useful metric for comparing the performance of different passive or noise-limited active 
sonars. For these cases, FOM is constant and independent of both range and environmental 
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propagation characteristics. It therefore captures all of the sonar design-related parameters 
for a given sonar operating in a specified noise environment against a defined target. 

Physically, FOM represents the maximum allowable one-way propagation loss resulting in 
zero signal excess. Wagner et al. note that the FOM metric for active sonar “today has fallen 
into disuse in both shipboard detection range prediction [and] in sonar system design,” 
since the benefit of a constant metric is lost in the reverberation-limited case when FOM 
must vary with range and environment. “Instead, the individual components of active sonar 
equations ... are each considered separately in arriving at a performance prediction” [84]. 
Since analysis of this type often requires a sophisticated MOS like CASS, however, we 
have elected to continue using FOM and limit our analysis to the noise-limited case. 

Assuming target detection is possible when SE > 0, we compute FOM by substituting 
Equation (2.2) into Equation (2.3) when SE = 0. Combining terms yields an expression 
for FOM as a function of the relevant design parameters [54]: 

0 = SL- 2EOM + TS-iN-DI+10 logjo B) - (5 logjo d - 10 log^o BT - 5 logjo n) 
FOM = ^(SL + TS-N + DI+10 logjo T - 5 logjo J + 5 logio • (2-4) 

We briefly describe each of these parameters, and provide sample calculations for the val¬ 
ues used in our simulations. Parameter values specific to the models corresponding to 
individual sonar designs are derived in Section 2.6.1, Section 2.6.2, and Section 2.7. 

• Directivity index (DI) of a transducer array describes its ability to “concentrate trans¬ 
mitted sound in a given direction” (Dlt), and improve the signal-to-noise ratio (SNR) 
received from a given direction (Dir). The DI is “a special case of the array 
gain (AG), where the signal is coherent and the noise is incoherent. This parame¬ 
ter is simpler to calculate and will normally be a satisfactory measure of the increase 
in SNR due to the array” [54]. This parameter is a function of the sonar’s design 
frequency and array geometry. 

• Source level (SE) of a projector array is a function of its acoustic power, P. If the 
array is directional, SE also depends upon the Dlf. 
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• Pulse duration (T) determines a sonar’s range resolution, with shorter pulses pro¬ 
viding better resolution due to smaller eeho separation. A CW “pulse of eonstant 
frequency and duration T seconds” will have a bandwidth, B = IjT Hz. For FM 
pulses, “the frequency of the pulse changes during the T seconds duration of the 
pulse [and] the bandwidth, B, is not now the inverse of the pulse length” [54]. 

• Detection index (d)h used to determine a sonar receiving system’s detection thresh¬ 
old (DT), expressed as the SNR corresponding to preset values for probability of de¬ 
tection {Pj) and probability of false alarm (Pfa)- Here, Pd is the probability that am¬ 
plitudes measured at the receiver which exceed DT consist of signal plus noise; and 
Pfa is the probability they consist of noise only. This relationship is typically plotted 
as a function of 5 logjo d ona curve of receiver operating characteristics (ROC). The 
sonar models derived below assume 5 log^g d = 10 dB, corresponding to Pd = 0.5 
and Pfa = 10“^ as given by Table 7.8 in [54] and used in the sonar design examples 
included in this reference. 

• The number of pings (n) contributing to a detection decision effectively reduces a 
sonar’s DT as more information is considered. The sonar models derived below as¬ 
sume 5 logjQ n = 3 dB, equivalent to processing four pings per decision; this value is 
used for all sonar design examples in [54]. 

• Target strength (TS) is the ratio of the intensity of the sound wave reflected by an 

underwater target, relative to the intensity of the incident sound wave from an ac¬ 
tive sonar pulse. This quantity, expressed in dB, is a function of sonar frequency, 
target size, geometry, and the angle of incidence between the sonar pulse and tar¬ 
get. The goal of our MCM search problem is to detect small mines on the sea 
floor with geometry approximated by finite cylinders of radius 0.1 meter and length 
1.0 meter, with hemispherical ends. For incident angles normal to the mine, TS 
is computed for a cylinder with radius a, length L, and wavelength A using the 
formula TS = logjQ The worst case TS, however, for incident angles strik¬ 

ing one of the two ends, is computed for a sphere with radius a using the for¬ 
mula TS = log^o (^) = -26 dB. Therefore, we use the conservative value of 
TS = -30 dB [54]. 
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• Noise (N) refers to the intensity level of the background noise that a desired signal 
must overcome to be detectable at the sonar receiver. Ambient noise in the ocean is 
due to three main sources: 1) noise in the sonar receiving system induced by thermal 
agitation of water molecules against the face of the hydrophone; 2) ambient sea noise 
from wind and waves, rain, distant shipping, and biological activity; and 3) vessel 
self-noise due to machinery, propellers, electrical interference, and fluid flow [54]. 
All of these sources contribute to the ambient noise level in different frequency bands, 
as illustrated by Figure ??. The average spectral level for thermal noise is depicted by 
the straight line in the bottom right corner of the plot. This line dominates other noise 
sources at the high frequencies used for mine hunting sonar, i.e., above 100 kHz, and 
can be computed for a given sonar frequency / Hz by the equation [53]: 


Ntherm = -75 + 20 logjQ / dB. (2.5) 

Even though exact design parameters for Navy sonars are difficult to obtain (and potentially 
classified) [84], a figure of merit (FOM) suitable for relative performance analysis can still 
be estimated from sonar design reference manuals [54] or commercial sonar specifications 
[91]-[93]. Furthermore, once the FOM for a given sonar problem is known, it is easy to 
compute the signal excess along a moving vehicle’s trajectory x(t), since it depends only on 
the propagation loss due to distance between vehicle and target. Consequently, calculating 
detection probability based on signal excess is especially attractive for sensor-based motion 
planning algorithms. 

2.3 Propagation Loss 

An acoustic pulse loses intensity as it propagates through the water, as the radiated power 
spreads throughout a larger and larger area. This can be modeled as spherical spreading, 
in which power is radiated in all directions, or as cylindrical spreading, in which power 
spreads outward between two planar boundaries. The pulse is also attenuated by absorption 
losses due to fluid viscosity and molecular relaxation of dissolved salts in seawater. Both 
spreading and absorption are functions of distance. For our signal excess calculations, we 
assume that propagation losses are due only to spherical spreading and absorption, “a use¬ 
ful working rule for initial design and performance comparisons” [54]. In fact, this same 
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Figure 2.6. Average ambient noise spectra. Source: [90], from [53]. 
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assumption is often required “in order to make reverberation models feasible” [56]. Ne¬ 
glecting other loss sources (e.g., scattering and refraction) via the simplifying assumptions 
listed in Section 2.1, propagation loss is 

PL(x(t),co) = 201ogio (||d> - f(Oll) +a\\(o- f(0||, (2.6) 

where a is the frequency-dependent attenuation coefficient of seawater. Tabulated values 
of a can be found in sonar design references such as [54] and [55]. While a varies with 
depth, salinity, and temperature, it depends most strongly on the sound frequency. We 
therefore compute this parameter using an equation in [53], which estimates a as a function 
of frequency in kHz: 

0 11 44 f 2 

a = ' \ + - - —^ + 0.0003/2 ^ Q QQ 3 dB/km. (2.7) 

2.4 Instantaneous Detection Rate 

In search theory, “the detection rate approach to computation of detection probabilities has 
proved to be more robust than the geometric models” used by “cookie cutter” sensor mod¬ 
els [51]. Originally developed in [36], this method assumes that a sensor has a detection 
rate y(t) called the ''instantaneous probability density (of detection).” This rate may vary 
with time due to the motion of searchers and targets, or to reflect changing environmental 
conditions, for example. Continuously searching over a small time interval At constitutes a 
single glimpse or scan with the sensor. Each glimpse provides a detection opportunity with 
the instantaneous probability of detecting a target given by y(t)At. This leads to the well- 
known exponential detection model, discussed in Section 3.4, which quantifies detection 
probability as a function of time. 

In order to use the exponential detection model, we must first compute detection rates for 
our sonar models. Detection rates based on our noise-limited signal excess model vary with 
distance between a target location m and a search vehicle following the trajectory x{t). If 
we also assume that the signal excess in Equation (2.3) is a normally-distributed random 
variable with mean SE and variance the instantaneous probability of detection for a 


35 



single glimpse with a sonar ean be written using its eumulative normal distribution O [51]: 


= O 


SE 


cr 


( 2 . 8 ) 


Based on our seleetion of the detection index d in Section 2.2, the instantaneous detection 
probability p{t) = 0.5 when SE{t) = 0, meaning the sonar has an equal probability of 
detecting or missing a mine. Regarding the selection of cr, Washburn notes that “most 
practitioners use a value of cr somewhere between 3 and 9 dB for sonar detection in the 
ocean” [51]. A value of cr = 5.6, computed by adding typical variance values for each 
term in the sonar equation, is provided in [84]. Moreover, a study which used the Navy’s 
CASS/GRAB software to simulate acoustic mine detections under varying environmental 
conditions, including seasonal sound speed profiles, unknown wind speeds, and imprecise 
bottom types, observed signal excess variations of 3 dB, 6 dB, or 9 dB in most cases [57], 
[58]. Figure 2.7 plots the instantaneous probability of detection vs. signal excess for these 
three values of cr. 


Detection Probability vs. Signal Excess 



Signal Excess [dB] 

Figure 2.7. Instantaneous detection probability vs. signal excess and cr. 


36 

























































































To compute a sonar’s detection rate y{t) from its instantaneous detection probability p{t), 
we further assume that detection opportunities (glimpses) can be modeled as a Poisson pro¬ 
cess and occur with mean rate A. The so-called Poisson Scan model produces the detection 
xdXcy{x{t),i6) = Ap(x{t),io) [51]. 

2.5 Detection Performance Modifiers 

The sonar detection model developed thus far, based on signal excess remaining after sub¬ 
tracting propagation losses from a given sonar’s figure of merit, is omni-directional. It 
depends only on the distance between the sonar and a MLO. This function could be used to 
construct a lateral range curve and corresponding sweep width for use in standard coverage 
planning algorithms [23]. Most actual sonar systems, however, are designed to transmit 
and receive with a specific beam pattern (see, e.g.. Figure 2.1 and Figure 2.2) and do not 
perform equally well in all directions. Actual detection performance depends not only on 
a sonar’s distance from a target, but also whether (and how long) it is pointed in the proper 
direction, at the proper place and time, to ensonify the target. The vehicle platform must 
maneuver to accomplish this. Conversely, high-resolution imaging sonars which rely on 
platform motion to methodically scan the sea floor (e.g., sidescan sonar) or construct long 
virtual hydrophone arrays (e.g., synthetic aperture sonar), require precise navigation along 
straight line trajectories to generate accurate imagery [13], [94]. Excessive platform motion 
can often yield poor performance for these systems. It is clear that sonar performance is 
tightly coupled to a vehicle’s motion. In the following sections, we derive simple shaping 
functions to model a sonar’s three-dimensional beam geometry as well as its dependence 
on platform motion. 

2.5.1 Field of View Considerations 

To more accurately estimate a sonar’s true detection performance when mounted on a vehi¬ 
cle, we must enforce its actual beam geometry in three dimensions. First, we define angular 
limits for the sonar’s horizontal and vertical fields of view (FOV). These FOV boundaries 
exist in the sonar reference frame, but we will assume without loss of generality that this 
frame is identical to the vehicle’s body-fixed reference frame. Next, we calculate the vec¬ 
tors between the sonar and each potential target, and resolve them in the body-fixed ref¬ 
erence frame to compute the angle to each target relative to the sonar’s FOV. Finally, we 
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apply a shaping function to degrade deteetion performanee for targets that fall outside of 
the angular FOV limits. 


The horizontal plane geometry at an instant in time is shown in Figure 2.8 for a SeaFox 
USV and forward-looking sonar with horizontal FOV, apov, of 200 degrees. Parameters 
defined or resolved in the body-fixed frame are denoted with a superseript b. The positions 
of the USV and potential mine target loeations are defined in the inertial 

referenee frame {n}. We eompute the lower and upper limits on azimuth angle for a sonar’s 
horizontal FOV by the expressions 


b 


aL = - 


O^FOV 

2 


, and 


b O'fov 

O^u = 


(2.9) 


The veetor between the sonar and a target of interest in the inertial frame is then defined 
di% h.x = \ (Ox - X, ojy - ij\ = {dx, di)\ .To determine the azimuth angle to this target in 
the sonar’s FOV, the veetor Ax must be resolved in the body-fixed referenee frame using 
the vehiele’s heading angle if/, producing the body-fixed eomponents: 

^dx= "dxcos(if/) + "dy sin(if/), and (2.10) 

^dy = - "dxsm(if/) + "dycosiif/). (2.11) 

Then, using the four-quadrant inverse tangent, we eompute 

= atan2(^dy, ^dx). (2.12) 

In the same manner, we eompute the lower and upper limits on elevation angle for the 
sonar’s vertieal FOV as: 


^Sl = SdE - 
^Su = Sde + 


^FOV 

2 

^FOV 


, and 


2 


(2.13) 


Here, is a fixed downward elevation angle seleeted to ensure that the sonar ean en- 
sonify the sea floor. Some sonar systems are eapable of eleetronieally steering their beams 
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Figure 2.8. Instantaneous geometry for a forward-looking sonar’s horizontal 
field of view. 


to a specified but this angle is frequently determined by a fixed mechanical mounting 
angle. For a vehicle traveling in the horizontal plane at constant altitude h above the bot¬ 
tom, the elevation angle between the sonar and a mine on the sea floor is identical in both 
reference frames, computed by 


e = 


s = arctan 


-h 




arctan 


-h 


( 


(2.14) 


We now define scalar shaping functions which degrade sonar detection performance for 
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mines outside the sonar’s horizontal (or vertieal) FOV. Eaeh shaping funetion is eon- 
strueted from two logistic functions [95]. These S-shaped sigmoidal curves [96] are cali¬ 
brated to smoothly transition a scalar value from 1 to 0 at the angular limits of the sonar’s 
FOV. This value modifies the probability of detecting a mine at a given location based 
on signal excess only, thereby preventing detection if the mine location lies beyond these 
angular limits. The azimuth and elevation shaping functions are 


FLS 


Fa{x{t),Oj) = 


1 


-I- 


1 


1, and 


-I- gPa{‘’aL - *0") + ^Pa{‘’a - ^au) 


(2.15) 

(2.16) 


The Pa and parameters can be tuned to adjust the slope of their respective sigmoidal 
curves, as discussed in Section 2.5.3. Figure 2.9 plots the azimuth scale factor vs. target 
azimuth angle and several values of pa for an FLS with a 120-degree horizontal FOV. 


Azimuth Scale Factor 



Figure 2.9. Fa vs. azimuth angle and Pa for a forward-looking sonar with 
120-degree horizontal field of view. 
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Similarly, Figure 2.10 plots the elevation scale factor vs. target elevation angle and several 
values of for a sonar with a 30-degree vertical FOV mounted at sde = -15 degrees. 
Although the x-axes of Figure 2.9 and Figure 2.10 are labeled in degrees, both functions 
are actually computed for azimuth (or elevation) angles given in radians. 


Elevation Scale Factor 



Figure 2.10. vs. elevation angle and ps for a forward-looking sonar with 
30-degree vertical field of view mounted at -15 degrees. 


The plots in Figure 2.11 and Figure 2.12 illustrate how detection probability based on 
signal-excess can be shaped using a sonar’s FOV. The color map indicates the probability 
of detecting a mine relative to the sonar’s reference frame. Figure 2.11 shows the horizontal 
plane geometry. The left plot shows the omni-directional detection probability, based on 
signal excess, for a forward-looking sonar with FOM = 72 dB and cr = 9 dB. The right 
plot shows the modified detection probability after applying an azimuth scale factor F^ 
corresponding to a horizontal FOV of 120 degrees with pa = 10. Similarly, Figure 2.12 
shows the vertical plane geometry. The left plot shows the modified detection probability 
after applying an elevation scale factor F^ corresponding to a vertical FOV of 30 degrees 
with Pe = 40. The sonar head is mounted on a surface craft with sde = -15 degrees. Note 
that the different scales used for distance (x-axis) and depth (z-axis) distort the apparent 
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Figure 2.11. Omni-directional detection probability before (left) and after 
modification by Fa (right) for a 120-degree horizontal field of view. 



N orth [100m] N orth [100m] 

Figure 2.12. Detection probability after modification by for a 30-degree 
vertical field of view mounted at -15 degrees (left) and close up (right). 


beam angle in this plot, but a close up version shown at right, plotted with equal axis 
scaling, reflects the expected vertical FOV. 

For a sidescan sonar comprised of dedicated port and starboard arrays, it is still possible to 
construct a continuous shaping function that describes both fields of view over the entire 
range of azimuth angles 6 [-;7r,;r]. In this case, we define the lower and upper azimuth 
limits relative to the center of the starboard array’s FOV, amid = 7r/2 as 

b O'FOV , b , O^FOV 

~ ^mid 2 3.nd Off/ — Q^fjiid "I" 2 ‘ (2*17) 
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The shaping function for the starboard array is calculated as before, i.e., Fa(x(t),6j) = 

Fa(x(t),oj), using Equation (2.15) and these new angular limits. Next, the shaping 
function for the port side array is calculated as 




+ 


Combining these shaping functions yields 


(2.18) 


^^^FM{t),a>) = ^'^‘^FM{t),c^) - P^'-^FM{t),c^). (2.19) 


Figure 2.13 plots the azimuth scale factor vs. target azimuth angle and several values of pa 
for a SSS with a 10-degree horizontal FOV. 


Azimuth Scale Factor (Sidescan Sonar) 



-180 -135 -90 -45 0 45 90 135 180 

Azimuth Angle to Target [deg] 


Figure 2.13. Fa vs. azimuth angle and pa for a sidescan sonar with 10-degree 
horizontal field of view. 


Using the Poisson Scan model, we can now compute the instantaneous detection rate for 
a given sonar by combining the effects of signal excess, three-dimensional FOV geometry, 
and average scan rate in the expression 


y(x{t),cd) = Ap(x{t),(5)Fa(x(t),aj)Ff;{x(t),co). (2.20) 
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2.5.2 l\irn Rate Considerations 

As already mentioned, some types of sonar require stable, straight line motion of a vehicle 
platform to produce high-resolution imagery. Sidescan sonar (SSS), for example, stacks the 
backscattered signals received from successive pings to produce an image of the sea floor. 
The across-track dimension of the resulting image corresponds to the two-way travel time 
of each ping, while the along-track dimension is formed by the vehicle’s forward motion. 
Turning maneuvers, therefore, have a direct impact on SSS performance [97]. In fact, 
“yawing motions ... are considered to have potentially the most serious degrading effects 
on sidescan images, because yaw causes the beam footprint to move along-track a distance 
proportional to the distance across-track” [98]. We model this behavior by applying another 
scale factor to degrade detection probability as a function of the vehicle turn rate, r(t). We 
select the Gaussian-like expression 

( 2 . 21 ) 

This function reaches a maximum value of one for straight line motion, e.g., when r{t) = 0, 
but falls off smoothly for non-zero turn rates. The slope of this curve can be adjusted via 
the tuning parameter cr^, as shown in Figure 2.14. Using this scale factor, we compute the 
instantaneous detection rate for a sidescan sonar by the expression 

yix{t),6j) = Ap(x(t),oj)Fa(x(t),6j)FE(x(t),oj)Fr(x(t)). (2.22) 
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Turn Rate Scale Factor 



Turn Rate [deg/s] 

Figure 2.14. iv vs. turn rate r(t) and cr^. 


2.5.3 Numeric Considerations 

A primary consideration when seleeting the shaping funetions deseribed in Seetion 2.5.1 
and Seetion 2.5.2 is their numerieal smoothness. We shall see in Chapter 3 how instan¬ 
taneous deteetion rate ean be used to ereate an objeetive funetion for our optimal seareh 
problem. Having a smooth (i.e., differentiable) objeetive funetion is extremely helpful 
when performing numerie optimization. Another eonsideration is the ability to derive and 
eneode analytie expressions for the objeetive funetion gradients. The SNOPT software 
paekage used to solve our optimal seareh problem, for example, “is able to estimate gra¬ 
dients by finite differenees ... for eaeh variable whose partial derivatives need to be esti¬ 
mated. However, this reduees the reliability of the optimization algorithms, and it ean be 
very [eomputationally] expensive if there are many sueh variables” [99]. 

These shaping funetions were also designed to be flexible, as the logistie funetions ean 
be ealibrated to refleet most sonar FOV geometries by setting the azimuth and elevation 
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angular limits and tuning a single growth rate parameter: pa for azimuth angle, or for 
elevation angle. Eaeh parameter multiplies the exponents in Equation (2.15) and Equa¬ 
tion (2.16), respeetively, to eontrol the slope of its sigmoidal eurve. This defines the “erisp- 
ness” of a sonar’s EOV boundary between regions of high and low deteetion probability. 
Eow parameter values result in a gradual transition. As values inerease, however, this tran¬ 
sition tends toward a diseontinuous step funetion, whieh presents numerieal diffieulties. 
Therefore, we have derived a heuristie to guide the seleetion of an appropriate growth rate 
value based on two qualitative metries: 

• Scale Factor Threshold (SFT): the value the seale faetor should attain within the 
sonar’s field of view (EOV). 

• Fraction Below Threshold (FBT): the portion of the nominal EOV that is below the 
desired SFT. 


The logistie funetion, evaluated at a EOV boundary, ean be rearranged to ealeulate growth 
rates that satisfy these metries for the azimuth and elevation shaping funetions via the 
expressions 


21n(l - SFT) 
FBT- Hfov 


and Pe = - 


21n(l - SFT) 
FBT-Vfov ’ 


respeetively. 


(2.23) 


One ean experiment with these metries to arrive at a growth rate value that aehieves a 
balanee between realistie EOV boundaries and a smooth objeetive funetion for numerie 
optimization. 


Einally, we aeknowledge that for distanees less than one meter, the spherieal spreading 
term 201ogiQ (||m - T(0||) in Equation (2.6) will eontribute a negative propagation loss, 
sinee spreading losses are defined relative to an intensity measured one meter away from 
the souree. Even more eoneeming is the faet that this term is undefined when the distanee 
equals zero. However, sinee our MCM seareh problem is foeused on finding bottom mines, 
the distanee to any mine target is guaranteed to exeeed one meter as long as our seareh 
vehiele travels at an altitude of one meter or greater. 
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2.6 Forward-Looking Sonar Models 

This section derives figure of merit (FOM) values for the forward-looking sonar models 
used in this dissertation. We consider two different designs: 

1. A long-range, low-resolution sonar designed with a cylindrical transducer array to 
provide a wide horizontal FOV. This type of sonar is typically used for wide-area 
surveys to detect MLOs during the first phase of an MCM operation. 

2. High-resolution, blazed array imaging sonar suitable for reacquisition and identifica¬ 
tion of previously detected targets during follow-on MCM missions. 

Both FLS designs are examples of “sectorscan sonar,” which generate two-dimensional 
images from each pulse [94]. Examples of this imagery are shown in Figure 2.15. 



Figure 2.15. Sample images from a cylindrical array (left) and blazed array 
(right) FLS. Image at left is courtesy of Thunder Bay 2010 Expedition, 
NOAA-OER. Source: [100]. 


2.6.1 Cylindrical Array Model 

A cylindrical array of transducer elements is a common, practical sonar design. Individual 
elements are grouped in vertical lines, or staves, to obtain a desired vertical beamwidth, and 
multiple staves are arranged in a ring to provide the required azimuth coverage [54]. Arrays 
of this type can be found on submarines, as pictured in [58] and [101], and in systems like 
the Autonomous Topographic Large Area Sonar (ATLAS), shown mounted on the NFS 
SeaLox USV in Ligure 2.16. A LOM for a long-range detection sonar similar to ATLAS 
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Figure 2.16. An ATLAS forward-looking sonar mounted on the NPS SeaFox 
USV. 

can be computed using the design example for a mine-hunting sonar described in Chapter 
11 of [54], 

We specify a 200 kHz sonar with 120-degree horizontal field of view, 5-degree vertical 
field of view, and a nominal operating range of 400 meters. We assume this sonar transmits 
an FM pulse with bandwidth 5 = 80 kHz and duration T = 10 ms, which yields better 
noise-limited performance than a CW pulse for the values used in [54]. We further assume 
that the sonar’s projector array is comprised of multiple transducer elements so it can steer 
its beam in the vertical plane. From the expression in Table 2.3 of [54], the number of 
elements required for a sonar stave to achieve a beamwidth of 5 degrees is n = 100/BW = 
100/5 = 20. Assuming the horizontal transmit beamwidth is 120 degrees, the projector 
array requires only m = 1 vertical stave, and the transmit directivity index for this baffled, 
cylindrical array can be calculated from Table 2.5 in [54] by the expression 

DA = 3 + lOlogiomn = 3 + lOlogjofl • 20) = 16 dB. (2.24) 
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If we assume that the total acoustic power radiated by this projector is P = 10 Watts, we 
can compute the source level (SL) for the sonar as described in [54] by the expression 

SL = 10logio P + 170.8 + DIt = 197 dB. (2.25) 

Turning attention to the sonar’s receive array, we specify narrow, 2-degree horizontal and 
vertical beam widths so the sonar can resolve small MLOs in its field of view. Table 
2.5 in [54] provides a formula for calculating the receive directivity index of a baffled, 
cylindrical array based on its height h in meters, diameter d in meters, and design frequency 
/ in kHz. For this sonar, the receive directivity index is 

Dir = lOXogiQShdf. (2.26) 

Assuming half-wavelength spacing of its transducer elements (a function of the design 
frequency), the array’s height is h = 16/{BW^ ■ /) = 76/(2 • 200) = 0.19 m, while its 
diameter ish = SS/(BWh • /) = 88/(2 • 200) = 0.22 m [54]. Substituting these values into 
Equation (2.26), we compute Dir = 39 dB for this receive array. 

Next, we compute the attenuation coefficient using Equation (2.7), and the noise back¬ 
ground due to thermal agitation using Equation (2.5), both functions of the sonar’s 200 kHz 
design frequency. The attenuation coefficient is a = 52 dB/km and the noise due to thermal 
agitation is Ntherm = 31 dB. However, we use a more conservative value of A = 34 dB 
to compute figure of merit in Equation (2.4). This 3 dB increase can accommodate addi¬ 
tional self-noise from the vehicle platform at levels comparable to the calculated Ntherm 
value [54], e.g., for loud vehicles like the NFS SeaEox USV. 

Einally, we estimate the Poisson Scan rate for this sonar model using the concept of “ping- 
to-ping overlap” illustrated in Eigure 2.17. This capability generates multiple looks, from 
different viewpoints, at MEOs on the sea floor, providing better detection and mapping per¬ 
formance [100]. It also makes sectorscan sonar more robust to turning maneuvers than SSS, 
provided the vehicle platform has sufficiently accurate navigation. To ensure 95% ping-to¬ 
ping overlap from a sonar with 400-meter swath width, mounted on an AUV traveling 1.5 
m/s, the sonar must ping about every 10 seconds. Therefore, To accommodate faster USV 
platforms, we assume that our sonar model pings every 5 seconds, which corresponds to a 
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Figure 2.17. Ping-to-ping overlap for a wide-sector FLS. Image courtesy of 
Thunder Bay 2010 Expedition, NOAA-OER. Source: [100]. 

Poisson Scan rate A = 0.2 scans per second. 

The parameters used to model this sonar are listed in Table 2.1. All specified (or assumed) 
values are italicicized, while calculated values are listed in plain text. The operating fre¬ 
quency, projector source level, cylindrical array geometry, and pulse characteristics chosen 
for this design yield a FOM of 72 dB. 

2.6.2 Blazed Array Models 

A relatively recent sonar design technique, based on “blazed” acoustic arrays, has led to 
a class of smaller, lighter, lower power imaging sonars which are well-suited for deploy¬ 
ment from small autonomous underwater vehicles. Leveraging techniques from the fields 
of radar and optics, a blazed array can “map angular image information into the frequency 
domain” [102]. In principle, these acoustic arrays are analogous to optical diffraction grat¬ 
ings which can separate a broad spectrum signal (white light) into individual, angularly- 
separated frequencies (colors) [103]. Blazed sonar arrays separate a broadband acoustic 
pulse into a “frequency-dispersed sound field” in which each frequency corresponds to a 
separate sonar beam. Unlike traditional sonar designs that use dedicated electronics to 
form and steer the beams generated by each stave in the aray, e.g., [104], “this approach 
allows multiple independent beams to be simultaneously formed from a single hardware 
channel” [102]. 
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Table 2.1. Sonar design parameters used to calculate noise-limited FOM. 




Forward-Looking Sonar 

Sidescan 


Sonar Design Parameters 

Cylindrical Array 

Blazed Array 

Sonar 


Frequency 

200 kHz 

450 kHz 

900 kHz 

900 kHz 


NominalRange 

400 m 

200 m 

100 m 

40 m 


Transmit Beam Widths 

Horizontal 

120° 

o 

0 

o 

o 

o 

T- 

0 


Vertical 

5° 

10° 

20° 

0 

O 

s 

s 

Receive Beam Widths 

Horizontal 

2° 

1° 

1° 

o 

T- 

o 

o 

s 

Vertical 

2° 

10° 

20° 

O 

O 

’xl- 

Pulse Length, T 

FM 

10 ms 

10 ms 

10 ms 



CW 

— 

— 

— 

6.67 jxs 


Pulse Bandwidth, B 

FM 

80 kHz 

80 kHz 

80 kHz 



CW 

— 

— 

— 

150 kHz 


Detection Index, SlogjQ J 

lOdB 

lOdB 

lOdB 

lOdB 


Detection Pings, 5 logjQ n 

3 dB 

3 dB 

3 dB 

3 dB 


Poisson Scan Rate, A 

0.2 scan/s 

0.5 scan/s 

1.0 scan/s 

25 scan/s 


Attenuation Coefficient, a 

52 dB/km 

104 dB/km 

287 dB/km 

287 dB 


Directivity Index, DI 

Transmit 

16 dB 

24 dB 

24 dB 

27 dB 

c3 

Receive 

39 dB 

24 dB 

24 dB 

27 dB 

o 

Source Level, SL 

197 dB 

207 dB 

206 dB 

204 dB 

U 

Ambient or Self-Noise, N 

34 dB 

41 dB 

45 dB 

44 dB 


Target Strength, TS 

-30 dB 

-30 dB 

-30 dB 

-30 dB 


Figure of Merit, FOM 

72 dB 

66 dB 

64 dB 

49 dB 


In this section, we apply the procedure described in Section 2.6.1 to compute FOM values 
for two blazed array multibeam imaging sonars. Teledyne BlueView’s P450 Series and 
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Figure 2.18. Combined beam pattern for a 50-degree horizontal field of view 
FLS constructed from two individual blazed arrays. Source: [105]. 


P900 Series systems are modular designs eomprised of multiple blazed arrays that operate 
at 450 kHz and 900 kHz eenter frequencies, respectively. Figure 2.18 illustrates how in¬ 
dividual staves can be combined to form a larger FOV. We will model the P450-90 and 
P900-90 sonars which utilize four staves to produce a 90-degree horizontal FOV. 

An individual stave produces a 25-degree fan of beams in the image-plane, each with a dis¬ 
tinct frequency and angle relative to the face of the stave. The lowest frequency beam angle 
is 45 degrees, and the highest frequency beam angle is 70 degrees, as shown in Figure 2.18 
(left). The combined FOV for a two-stave system is depicted as a three-dimensional solid in 
Figure 2.18 (right), illustrating how the beam widths vary with frequency in both the image- 
and cross-image planes. Although beam pattern geometries for multi-stave systems like 
those depicted in [105] are complex, we assume that image processing algorithms allow us 
to model these sonars as conventional line arrays operating at the center frequency of their 
broadband pulse. Under this assumption, we can use manufacturer specifications for oper¬ 
ating frequency, field of view, beam width (in the image- and cross-image planes), number 
of beams, and update rate to calculate a FOM for the P450-90 and P900-90 sonars [91]. 
We further assume that the these sonars use an FM pulse duration T = 10 ms. 
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Note that the manufaeturer’s speeifieations list the total number of beams in a given sonar. 
Since these sonars are constructed using modular staves, each with 128 individual beams, 
we compute the directivity index for a single stave using the expression for a baffled line 
array with n = 128 elements from Table 2.5 in [54]. For this sonar, the directivity index is 

Dl, = D/, = 3 + lOlogion = 3 + 101ogio(128) = 24 dB. (2.27) 

The transmit source levels for the P450-130 and P900-45 sonars are given in [106] as 
207 dB and 206 dB, respectively. Using Equation (2.25) and the calculated value for 
DIt = 24 dB, we can solve for a total acoustic power level P between 13 and 17 Watts, 
which is reasonable given that the stated electrical power consumption for these models is 
between 15 and 30 Watts. 

Next, we use Equation (2.7) to compute attenuation coefflcients for the 450 kHz and 
900 kHz operating frequencies of a = 104 dB/km and a = 287 dB/km, respectively. Erom 
Equation (2.5) we estimate that the noise due to thermal agitation is Ntherm = 38 dB and 
Ntherm = 44 dB, respectively. After accommodating additional self-noise noise due to the 
vehicle platform (which has lesser impact at higher frequencies) we use conservative values 
of A = 41 dB and A = 45 dB, respectively. 

Einally, we estimate the Poisson Scan rate for both sonar models by scaling the maximum 
update rate specified by the manufacturer. Assuming that these listed values apply to a 
sonar operating at its minimum optimal range, we scale the listed values by the nominal 
operating ranges used in our problem. Eor the P450-90 and P900-90 sonars, we estimate 
scan rates of T = 0.5 and T = 1.0 scans per second, respectively, based on nominal operat¬ 
ing ranges of 200 and 100 meters. These values agree with practical update rates observed 
when deploying these sensors on a REMUS AUV. The parameters used to model the P450- 
90 and P900-90 blazed array sonars are listed in Table 2.1 alongside the parameters for the 
cylindrical EES. The resulting EOM values calculated for these models are 66 dB and 
64 dB, respectively. 

2.7 Sidescan Sonar Model 

Next, we estimate a figure of merit for a short-range, side-looking sonar similar to the 
sidescan sonar used on the NPS REMUS 100 AUV. This type of sensor is representative 
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of high-resolution sonars used to reacquire previously-detected MLOs and identify them 
for subsequent neutralization [13]. Following the design example for a SSS in Chapter 
10 of [54], we use manufacturer specifications to derive a model for a 900 kHz sidescan 
sonar [92]. In contrast with the forward-looking sonars considered in Section 2.6, this sonar 
transmits a CW pulse with duration T = 6.67 yus and bandwidth 5 = 1/r = 150 kHz. It 
also has a very narrow horizontal FOV and a wide vertical FOV. As before, we can use 
the expression in Table 2.3 of [54] to compute the number of elements in a line array from 
its beamwidth n = 100/5VF = 100/0.4 = 250. Substituting n into Equation (2.27), we 
calculate the directivity indices for the SSS’s transmit/receive arrays: Dlf = = 27 dB. 

Assuming the sonar radiates 4 Watts of acoustic power, the source level for this sonar 
is SL = 204 dB, using Equation (2.25). The attenuation coefficient and thermal agita¬ 
tion noise are computed from the 900 kHz operating frequency as a = 287 dB/km and 
Ntherm = N = 44 dB, respectively, since underwater platforms have much lower self¬ 
noise than surface craft. Because the ping rate of sidescan sonar is usually determined by 
the vehicle platform’s speed and the sonar’s range setting [107], we estimate a Poisson 
Scan rate for this sonar based on the two-way travel time of sound for a nominal operating 
range of 30 meters: A = 2 ^( 3 q = 25 scans per second. The parameters used to model this 
900 kHz SSS are listed in Table 2.1 and yield a EOM of 49 dB. 

2.8 Model Verification and Validation 

The sonar models developed in Section 2.6 and Section 2.7 were tested in simulation to 
verify their detection performance. The primary objective was to confirm whether a given 
sensor model’s signal excess, EOV geometry, and Poisson Scan rate are accurately repre¬ 
sented. Eor example, Eigure 2.19 illustrates the sensor coverage obtained by a USV with 
200 kHz EES as it traveled through a square search area along the open loop trajectory 
shown. In this plot, the color map represents the probability that the searcher failed to 
detect a mine at each map location before the end of the mission. This probability of non¬ 
detection, PNoiTf), conditioned over the entire search area, is our measure of MCM risk. 
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Figure 2.19. Simulated search with USV and forward-looking sonar. 


Assuming that we have no prior information and that a mine ean be uniformly plaeed 
anywhere within the seareh area, the maximum P^oiTf) value (dark red) is a funetion 
of the seareh area size. In areas eovered by the USV’s sonar, PNoiTf) approaehes zero 
(blue), indieating a high probability of detecting a mine in those areas. We note the swath 
width produced by the sensor model is approximately equal to twice the nominal range, as 
expected. 

A similar trajectory is shown for an AUV with sidescan sonar in Figure 2.20. To accommo¬ 
date the slower AUV and shorter range sonar, the designated search area is much smaller 
than the FLS search area in Figure 2.19. 
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Figure 2.20. Simulated search with AUV and sidescan sonar. 


As a result, the maximum values for P^oiTf) are higher in a given area, a funetion of the 
uniform probability being spread over a smaller area. Nevertheless, the relative color scale 
indicates that our sidescan detection model is consistent with observed sidescan perfor¬ 
mance, namely it provides little to no coverage in the near-nadir region beneath the vehicle 
and performance is severely degraded when the vehicle turns. 

In this chapter, we have defined detection models for several different active sonars, based 
on noise-limited signal excess. These models provide the instantaneous probability of de¬ 
tecting a target as a function of sonar design parameters, the figure of merit (FOM), and 
the range to the target, which determines propagation losses due to acoustic spreading 
and absorption. We then defined shaping functions which enforce a sonar’s three dimen¬ 
sional geometry based on its horizontal and vertical FOV. After applying the azimuth and 
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elevation angle shaping funetions, the Poisson Sean model is used to eonvert the instan¬ 
taneous deteetion probability into a deteetion rate. Chapter 3 will formulate the optimal 
search problem as a GenOC problem, and utilize these sonar detection models to generate 
sensor-based motion plans to accomplish various MCM missions with different vehicles 
and sensors. 


57 



THIS PAGE INTENTIONALLY LEET BLANK 


58 



CHAPTER 3: 

Optimal Search Problem Formulation 


In this chapter, we cast different mine countermeasures (MCM) operations as optimal 
search problems whose solutions yield motion plans for a team of autonomous vehicles. 
We specifically consider two MCM missions: 1) an initial wide area survey with long- 
range, low-resolution sonar to detect and locate MLOs, and 2) a subsequent mission to 
revisit these locations with high-resolution sonar for positive target identification. We for¬ 
mulate an objective function which can incorporate different searcher, sensor, and target 
distribution models to solve different MCM search problems in the GenOC framework, 
demonstrating its flexibility. The benefits of this approach are two-fold: solutions not only 
specify trajectories that each vehicle should execute for a given sensor payload, but also es¬ 
tablish performance benchmarks for a given problem. The former can determine the most 
effective search pattern for a given vehicle or sensor, as discussed in Chapter 4. The latter 
provides a quantitative baseline for comparing different system configurations, which we 
explore in Chapter 5. 

A number of vehicles and sensors are capable of performing the search tasks described 
above, but mission planners must consider which combinations are most effective for a 
given MCM operation. Often, the available vehicle platform dictates which sensors can be 
utilized, while a sensor may dictate the type of trajectory a vehicle must follow. We assume 
in this work that each vehicle deploys only one type of sonar, but acknowledge that some 
developmental systems can carry multiple sophisticated sonars at once [13], [108]. An¬ 
other characteristic of search operations is that prior information (or the lack thereof) about 
potential target locations influences how the search is planned and executed. The GenOC 
framework includes all of these characteristics: multi-vehicle operations, sensor-based mo¬ 
tion planning, and prior information. It can be customized to explore a wide variety of 
MCM scenarios simply by swapping different models of vehicle/sensor performance and 
initial target distribution. 
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The optimal search problem defined in this chapter assumes that: 

• Targets are bottom mines with known target strength. 

• Sea floor is flat. 

• Water depth is constant. 

• Search effort is confined to a rectangular area. 

• Available mission time is fixed. 

The following sections describe the mathematical models and the objective function used 
to solve MCM search scenarios within the GenOC framework. 


3.1 Searcher Models 

The search vehicles selected for a given mission place constraints on the admissible so¬ 
lutions to an optimal control problem. Specifically, we define a mathematical model for 
a vehicle’s dynamic equations of motion. This model relates the vehicle’s state variables 
x{t) 6 and control inputs u{t) 6 R^“ through a set of ordinary differential equa¬ 
tions (ODEs) expressed in state-space form: x{t) = / {x(t),u(t)). This model places 
dynamic constraints on how the states may evolve with time. Similarly, there may be 
algebraic constraints placed on the control inputs due to physical actuator limitations: 
umin < m(0 < Umax for all t 6 [0,r/]. Additional constraints can be specified in or¬ 
der to bound the region of state space explored during optimization: xmin ^ -^(0 ^ ^max 
for all t 6 [0,r/-]. We note that for multi-vehicle problems with Nv searchers, the state 
and control vectors are simply expanded to include the states and controls of each vehicle, 
x(t) 6 R^"^^ and u(t) 6 R^"^“, respectively. The following sections describe models of 
two autonomous vehicles which are representative of naval platforms used for MCM. 

3.1.1 SeaFox Unmanned Surface Vessel 

The NFS SeaFox USVs are small, 5-meter rigid hull inflatable boats (RfflBs) originally 
designed for remote-controlled intelligence, surveillance, reconnaissance (ISR), force pro¬ 
tection, and maritime interdiction operations (MIO) conducted by the USN and U.S. Coast 
Guard (USCG) [109], [110]. CAVR has converted these vessels into fully autonomous 
surface craft in support of various research programs, including sonar-based path plan¬ 
ning for riverine navigation [111], [112] and precise speed control [113]. More recently. 
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CAVR modified its SeaFox Mk II USV to deploy an ATLAS minehunting FLS for MCM 
research [114]. 

To develop a model for these vehicles, we assume that USVs conduct MCM search 
missions at constant velocity, without aggressive maneuvers, and therefore exhibit sim¬ 
ple planar motion at the sea surface (i.e., pitch, roll, and heave motions are zero). If 
we further assume that sway motions are negligible (i.e., sideslip is zero) the equa¬ 
tions of motion can be adequately modeled by kinematics only. Using the state vector 
x{t) = [v(0, , the state-space equations of motion (EOM) are: 

i(0 = V cos(i/r(0) 
y(t) = V sin(i//(t)) 

<A(0 = r(t) 

r(t) = ^ (Ku(t) - r(t)). (3.1) 

The state variable pair y{t)] defines the vehicle’s position in meters along the north 
and east axes of the inertial reference frame; i/^(0 describes the vehicle’s heading angle in 
radians measured clockwise from the North axis; and r{t) is the vehicle’s turn rate in ra¬ 
dians per second. The vehicle travels with constant forward velocity V meters per second, 
measured along the body-fixed x-axis. Equation (3.1) implements a first-order approxi¬ 
mation to the well-known Nomoto model for ship steering equations, a simple transfer 
function between rudder angle u(0 = 5r(t) and turn rate r(t) that “is widely used for ship 
autopilot design due to its simplicity and accuracy” [115], [116]. The Nomoto gain con¬ 
stant K in inverse seconds, and time constant T in seconds can be identified from sea trial 
maneuvers as described in [117]-[119]. Table 3.1 lists the values of V, K, and T used in 
our SeaEox USV model. 

Table 3.1. Design parameters for unmanned surface vessel model. 


Design Parameter 

SeaEox Value 

Nomoto Gain Constant, K 

0.5 1/s 

Nomoto Time Constant, T 

5.0 s 

Velocity, V 

2.5 m/s 
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3.1.2 REMUS 100 Autonomous Underwater Vehicle 

The REMUS 100 AUV is a small, rapidly deployable unmanned underwater vehicle for 
collecting environmental data in the ocean [120]. Its modular design accommodates a 
number of different sensors for hydrographic survey missions, and its SSS system can 
make detailed maps of the ocean floor. One of the first AUVs adopted for naval MCM 
operations [82], [121], REMUS vehicles were used during Operation Iraqi Ereedom in 
2003 [122]. The REMUS family of vehicles includes two MCM variants in use by the 
Navy today: the MK 18 Mod 1 Swordfish, based on the 7.5-inch diameter REMUS 100, 
and the MK 18 Mod 2 Kingfish, based on the 12.75-inch diameter REMUS 600 [123]. 
CAVR operates three REMUS 100 AUVs in support of its research programs, and has been 
developing sensor-based navigation algorithms that utilize blazed array forward-looking 
sonar since 2004 [124]-[126]. 

Autonomous underwater vehicles can move in all three dimensions, and six degrees of 
freedom (DOE) are required to describe this motion completely. An example of a full 
6 -DOE model for simulating the nonlinear dynamics of a REMUS 100 is presented in 
[127]. In practice, however, these EOM are usually decoupled into separate, linearized 
equations in the horizontal and vertical planes so that designers can develop controllers for 
steering and diving, respectively. Eor our search problem, since AUVs typically conduct 
constant-velocity SSS surveys at a fixed altitude above the bottom, we consider only two- 
dimensional planar motion. Einally, as a matter of convenience when implementing multi¬ 
vehicle problems in software, we prefer a motion model with the same form as the SeaEox 
USV model in Section 3.1.1. This provides easier state vector indexing when AUVs and 
US Vs operate in a heterogeneous vehicle team. 

Under the same assumptions of zero pitch, roll, and heave motion, we can derive a Nomoto 
steering model for the REMUS 100 AUV from the linearized, decoupled, lateral steering 
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equations in Equation (118) of [127], repeated here: 

m-Yi, -Yf- 0 v(t) -Yy muQ-Yr 0 v(t) Yg^ 

-No hz-Nr 0 r(t) + -N, -Nr 0 r(t) = Ns, Sr(t) (3.2) 

0 0 ij [iA(oJ [ 0 -1 oj [iA(oJ [ 0 

where uq = steady-state surge velocity in the x-direction, 
v(t) = sway velocity in the y-direction, 
m = vehicle’s mass, 

fzz = vehicle’s yaw moment of inertia about the z-axis, 

Y = linear hydrodynamic coefficients producing sway forces, and 
N = linear hydrodynamic coefficients producing yaw moments. 

In general, control inputs and state variables (and their derivatives) produce nonlinear hy¬ 
drodynamic forces and moments. It is common practice, however, to approximate these 
effects by multiplying each contributing variable with a linearized hydrodynamic coeffi¬ 
cient. In Equation (3.2), Y and N denote coefficients that produce sway forces and yaw 
moments, respectively, while subscripts denote their corresponding control input or state 
variable. Assuming that sway velocity is zero (no sideslip), we rearrange Equation (3.2) as 

Izz-Nr oj [f(oj \-Nr oj [r(0 

0 1 (//(t) -1 0 

{Izz-Nr)r{t)-Nrr{t) 

(hz - Nr) r(t) 
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Manipulating the r{t) expression from Equation (3.3) into the form of our first order 
Nomoto steering model (see Equation (3.1)), we have 


{h,-N,)m 

f(t) 

f(t) 

where: T 

K 

u(t) 


Nrr(t) + Ns,Sr(t) 

r(t) + 


(h^-NrY ■ ihz-Nr) 

Nr - hz 
Nr ’ 

Nr' 

Sr(t). 


and 


(3.4) 


Substituting values from [127] for the yaw axis moment of inertia Yz and the hydrody- 
namie eoefiieients Nr, Nr, and (using the SeaEox sign eonvention), we ealeulate the 
parameters for our REMUS 100 model listed in Table 3.2. 

Table 3.2. Design parameters for autonomous underwater vehicle model. 


Design Parameter 

REMUS 100 Value 

Nomoto Gain Constant, K 

2.0 1/s 

Nomoto Time Constant, T 

1.0 s 

Veloeity, V 

1.5 m/s 


Sinee our MCM seenario eoneems the seareh for bottom mines, the altitude h must also be 
speeified for a given vehiele’s mission. Eor surfaee eraft, altitude h equals the water depth 
under our flat bottom, eonstant depth assumptions. 

3.2 Sensor Models 

Physies-based models for different types of aetive sonar used in MCM were developed 
in Chapter 2. These models ealeulate the instantaneous probability p that a given sonar 
ean deteet an eeho from a speeifie target against the expeeted ambient noise level. This 
quantity is a funetion of the sonar design, parameterized by its FOM, and the two-way 
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propagation losses PL between the sonar and target when a sean oeeurs. At these sean 
times, the instantaneous deteetion probability also depends on whether the target lies within 
the sonar’s three dimensional FOV, a function of sonar geometry and vehicle trajectory 
x{t). Three scalar shaping functions were designed to characterize these relationships: Fa, 
Fs, and iy. For a given sonar, this process occurs at an average rate A, producing the 
instantaneous detection rate y defined in Equation (2.22). Including the terms described 
above yields 

y(x(t),oj) = AFa(x(t),oj)Fs(x(t),oj)Fr(x(t))p(x(t),co) 

= AFaim,^)Feim,(^)Frim)<^ (xU),co) 

= AFa(m,co)F,(m,co)Fr(m)^ ^ FOM-PL{x (0,cu) j 

The expressions used to compute each term in Equation (3.5) are listed in Table 3.3. 


Table 3.3. Terms used to compute instantaneous detection rate. 


Symbol 

Definition 

Cross Reference 

A 

Poisson Scan Rate 

Section 2.4 

SE 

Signal Excess 

Equation (2.3) 

FOM 

Eigure of Merit 

Equation (2.4) 

PL 

Propagation Eoss 

Equation (2.6) 

P 

Probability of Detection 

Equation (2.8) 

Fa 

Azimuth Shaping Eunction 
forward-looking sonar 

sidescan sonar 

Equation (2.15) 
Equation (2.19) 

Fe 

Elevation Shaping Eunction 

Equation (2.16) 

Fr 

Turn Rate Shaping Eunction 

Equation (2.21) 


3.3 Target Models 

The GenOC framework was developed to address optimal control problems with parame¬ 
ter uncertainty [5], [10], [11]. Eor the MCM search problems considered in this disserta- 
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tion, the uncertain parameter is the location of a mine target on the sea floor. We assume, 
therefore, that target location is a stochastic parameter co distributed over a search area Q 
according to a known continuous PDF: 0(m) : R. That is, m 6 c R^" [5], and 

No, = 2. We can model a number of possible target distributions by specifying a different 
PDF, and this directly influences the solution of the optimal search problem. 

We can use different target distribution models to simulate various MCM missions. In the 
initial phase of an MCM operation, for example, wide-area surveys are conducted to detect 
and localize MLOs in a given search area. We formulate this task as an optimal search 
problem where no prior target data is available. We therefore model the target distribution 
with a joint uniform PDF, bounded by the search area coordinates. This PDF contains 
no exploitable information, i.e., there is an equal probability of finding a target anywhere 
within the search area. 


In subsequent phases of an MCM operation, additional sorties are conducted to reac¬ 
quire previously detected MLOs and identify mines from non-mine/mine-like bottom ob¬ 
jects (NOMBOs) using high-resolution sonar. We can formulate this RID task as another 
optimal search problem which not only incorporates different vehicle and sensor models, 
but also leverages MLO location data gathered during a prior survey. The accuracy of prior 
information is commensurate with the navigational performance of the survey vehicle, so 
we model this variation with a PDF. Any continuous PDF can be specified. Walton, how¬ 
ever, suggests the use of joint normalized beta distributions for modeling purposes, due to 
their easy algebraic manipulation and customization (via the a, /3 shape parameters), as 
well as their finite radius of effectiveness [79]. The PDF for the beta distribution, defined 
for X 6 [0,1], a > 0, and J3 > 0 [128] is 


(()ix;a,/3) 
where B(a,/3) 


1 


Bia,/3) 
(a-mfi 


x“-i(l 




1 )! 


(or + yS- 1)! 


(3.6) 
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3.4 Objective Function 

In this section we present the exponential detection model, first described in [36], which is 
commonly used to quantify search performance in continuous time. Based on a sensor’s 
instantaneous detection rate (see Section 2.4), this model provides a convenient objective 
function for optimal search problems, with recent examples provided in [5], [6], [11]. For 
our problem, we define residual MCM risk as the probability that a team of autonomous 
vehicles fails to detect the mines in a search area by the end of an MCM operation. This 
scalar quantity can be readily calculated for a given set of vehicle and sonar capabilities, 
and also reflects the time available for search. Therefore, we utilize MCM risk as the 
objective function for our optimal search problems; minimizing this quantity maximizes 
the mission’s probability of success. 

Given an instantaneous detection rate yit), derivation of the exponential detection model 
proceeds from two key assumptions [36]: 

1. The probability of detection in the short time interval [t,t + At] is y{t)At. 

2. Detection events in all such non-overlapping time intervals are independent. 

Washburn cautions that the independence assumption may not hold in all situations. For 
example, consecutive detection failures due to low signal excess could be caused by low tar¬ 
get strength or poor acoustic conditions. Empirically, however, these assumptions “provide 
good approximations in a wide variety of circumstances” [51]. Koopman acknowledges 
the importance of recognizing when this assumption is legitimate or not, and justifies its 
use beyond cases of random search: 

The assumption is in fact legitimate—and important—when applied to condi¬ 
tional probabilities of detection: probabilities calculated on the basis of postu¬ 
lated positions and motions of the target. [37] 

This is precisely the case described by our objective function for MCM risk in Equa¬ 
tion (3.11), which we now derive using Koopman’s “assumption of independence.” 

Eet p{t) be the probability of detection at time t. Then, by the complement, the probability 
of a detection failure is q{t) = 1 - p{t). Under our stated assumptions, this probability 
becomes q{t -t- At) = q(t)(l - y(t)At) at the end of the next scan interval, which can be 
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rearranged as the differenee equation: 


q{t + At) - q{t) 

- 

In the limit as At ^ 0, we obtain the differential equation 

q{t) = -q{t)y{t). 


(3.7) 


(3.8) 


whieh has the elosed form solution 


q(t) = (39^ 

and leads to the exponential deteetion model: p(t) = 1 - q(t) = I - e~ -fo 

Equation (3.9) represents the probability that a target was not deteeted by time t, so the 
residual MCM risk after eompleting an operation of mission duration is: 

(3.10) 


The objeetive of our optimal seareh problem is to minimize this risk. However, the in¬ 
stantaneous deteetion rate in Equation (3.5) depends on the vehiele trajeetory x(t) and the 
uneertain target loeation to, a random variable defined in Seetion 3.3. Consequently, Equa¬ 
tion (3.10) is itself a random variable, whieh we eannot minimize explieitly. Instead, we 
must minimize its expeeted value, eonditioned on the PDE of the target distribution. Henee, 
the objeetive funetion for a single vehiele beeomes 

J = E{qiTF)}= [ (3.11) 

Ja 

Eor missions with multiple searehers, we assume that instantaneous deteetion rates are 
additive, but vehiele- and sensor-speeifie. Eor example, the eombined deteetion rate for a 
team eomprised of N,j vehieles is 


N„ 


r(f(o,m) = 


k=\ 


FOM^ -PL ix(t),co) 


tr" 


(3.12) 
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The additive rate in Equation (3.12) makes the implicit assumption that multiple sonars 
do not acoustically interfere with one another. This assumption is rarely realistic, except 
when sonar systems have widely-separated design frequencies or when they are deployed 
in different locations, but this can be addressed during mission planning by imposing min¬ 
imum separation distances. Interestingly, the exponential detection model often produces 
multi-vehicle motion plans which resemble manually-separated vehicle trajectories. Since 
this model yields diminishing returns when multiple vehicles search the same location [51], 
it tends to encourage multi-vehicle solutions which distribute search effort to explore new 
regions of the search space. Under our assumptions, the expected residual MCM risk after 
a multi-vehicle operation is 

= f (3.13) 

Ja 

The objective functions in Equations (3.11) and (3.13) are differentiable, analytic expres¬ 
sions. Although somewhat tedious, it is possible to derive formulas for their gradients 
with respect to the state and control variables. This has benefits when using gradient-based 
numerical optimization algorithms. Encoding these formulas as user-defined functions sup¬ 
plied to the SNOPT optimization package, for example, significantly reduces the run time 
required to compute an optimal solution [99]. The objective function gradients for a single 
vehicle are derived in Appendix A. 

We observe that our objective functional 7, representing the expected residual MCM risk, 
has the same form as the running cost in Equation (1.2), the objective functional of the 
GenOC problem discussed in Chapter 1. This objective functional is repeated here for 
comparison with Equations (3.11) and (3.13): 

7 = y E (^x(Tf),coi^ + G R (x(T),u(T),T,aj) dr^ (()(aj)dio. (3.14) 

The end-point cost E i^x{Tf), in Equation (3.14) has been omitted from our objective 
functions. Meanwhile, the function G(-) of the running cost derives from the exponential 
detection model, i.e., G(-) = e~^’\ and R{-) = y(-)- The objective functional in Equa¬ 
tion (3.14) has been used to solve optimal search problems with multiple searchers and 
moving targets in cases where target motion can be conditionally-determined by uncertain 
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initial conditions [5], [7], [81]. Using similar objective functions in Equations (3.11) and 
(3.13) allows us to leverage the mathematical and computational framework previously de¬ 
veloped to handle this class of parameter-distributed, nonlinear optimal control problems. 
A contribution of this dissertation, therefore, is the development of new physics-based 
sonar models which allow the GenOC solution framework to address real-world MCM 
mission planning and analysis problems. 


3.5 Problem Scaling 

We need to solve these optimal search problems numerically, but the domains of our state 
variables x{t), control inputs u{t), uncertain parameters co, and objective function J all 
have different orders of magnitude. The search area, for example, may cover several square 
kilometers, while the objective function evaluates to a probability in the range [0,1]. It 
is important, therefore, to properly scale the problem before unleashing a numeric solver. 
This can be achieved by defining canonical units for distance, time, etc. and transforming 
the original problem’s variables into non-dimensional versions with similar domains [2]. 
Several examples which use variable scaling to numerically balance the equations of an 
optimal control problem are provided in [74]. 


For our search problems, the vehicle models from Section 3.1 can be scaled by canonical 
units of distance DU, time TU, and velocity VU = DU ITU to produce dimensionless 
variables designated by overbar notation: 


X 

7 

7 

t 

u 


X 


DU 

y 

DU 




l/TU 

t 

TU 


= {TU)r 


u. 


(3.15) 


Note that angular variables for heading ip and control input u (rudder angle) are already 
expressed in dimensionless units of radians. The chosen scaling must also be applied to 
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constant model parameters sueh as veloeity V, as well as the gain K and time T eonstants 
of the Nomoto steering model: 


K = = {TU)K (3.16) 

T 

T = —. 

TU 

Substituting these expressions into our original expressions for x and y yields their state 
spaee equations in nondimensional units: 


X = 

X = 


dx d(DUx) 

It ~ d{TUt) 
1 1 


DU dx 


= VUx 


vu"" VU 


TU dt 
V cos(if/) = V oos(i/r). 


(3.17) 


Similarly, we have 


y = V sin((/r). 


(3.18) 


Likewise, sealing by eanonieal units for ij/ and r yields the non-dimensional expressions 


dif/ dif/ \ dij/ 1 — 

~~dt ~ d(TUt) ~ fu~^ ~ fu"^ 
if/ = TUij/ = (TU)r = r, and 


dr d{llTU)r 1 dr 
^ ~~dt~ d{TUt) ~ fu^H 

f={TU^)f = {TU^)^{Ku-r) 


= (TU^)-^ 
TUT 


\tU TU I 




1 ^ 

-r 

rt/2 


(3.19) 


(3.20) 


Equations (3.17) through (3.19) eonfirm that our sealing has not ehanged the underlying 
dynamies of the problem. To ensure the objeetive funetion J is caleulated properly, physieal 
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units in the detection rate equation must also be scaled by the appropriate canonical units. 
The Poisson Scan rate A is scaled using TU to yield the non-dimensional form 


A = = {TU)A. (3.21) 

Recall from Equation (2.6) that our range-dependent propagation loss includes a spherical 
spreading term, 201ogio (||di - ic(0||), and an acoustic absorption term a ||m - x(t)\\. At 
each time t we can compute the distance D between a vehicle and target. Scaling this 
distance yields: 


D = \\oj- f|| 


D = ^(a)x- x)2 -t (a)y - «/)2 -t (ojz - z)^ 


D = -^(dx)^ -1- (dy)^ -1- (dz)^ 

(3.22) 

D = ^(DU'^)^ + (DUdy)^ + {DUlz)'^ 


1 -2 -2 -2 

D = DU yj dx + dy + dz 

D = DUD. 

(3.23) 


So D = DjDU, as expected. The attenuation coefficient a has units of dB/m and must be 
scaled by canonical distance DU (dB represents a ratio and is dimensionless already) as 


a = 


= {DU)a. 


l/DU 

Finally, the propagation loss calculated using these dimensionless quantities becomes 


(3.24) 


PL = 20logio (D) + a{D) (3.25) 

PL = 20logio [dud] + ^ (dud) 

PL = 20logio i^U) + 20 logio P) + « P), (3.26) 

which has an additional term due to the distance scale factor. The level curves in Figure 3.1 
verify that propagation loss calculated for non-dimensional ranges with Equation (3.26) 
and DU = 100 meters (shown at right) are equivalent to propagation loss calculated for 
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Level Curves vs. Range Level Curves vs. Range 




Figure 3.1. Propagation loss vs. physical range (left) and non-dimensional 
range (right) for 200 kHz FLS with figure of merit of 72 dB. 

physical ranges with Equation (3.25) (shown at left). These eurves also show the signal ex- 
eess vs. range for the 200 kHz FLS with FOM = 72 dB, derived in Seetion 2.6.1. Note that 
signal exeess is zero at a range of 400 meters. As an example, suppose we want to solve an 
MCM seareh problem for the SeaFox USV (see Seetion 3.1.1) eondueting a mine deteetion 
survey with the 200 kHz FLS above. Typieal bounds on the states, eontrols, seareh area, 
and eonstant parameters are defined in Table 3.4. This table ineludes their physieal val¬ 
ues (before sealing), and their non-dimensional values after sealing by the eanonieal units 
DU = 100 meters, TU = 100 seeonds, and VU = ^ = 1 m/s- 

3.6 Feasibility 

Most optimal eontrol problems eannot be solved analytieally. Often, numerieal methods 
are required to generate “optimal” trajeetories of the state variables and eontrol inputs that 
minimize a desired objeetive funetion, subjeet to eonstraints defined by the user. We must 
remember, however, that numerie solutions are ealeulated for a discretized version of the 
original problem. They meet the mathematieal definition of feasibility as long as all of 
the problem eonstraints are satisfied at a finite number of nodes eomprising the diserete 
problem [129]-[131]. 
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Table 3.4. Example of physical and non-dimensional parameter domains. 


Parameter 

Physical Domain 

Canonical Domain 

Time 

0 <t < 1800 

s 

0 < / < 18 

Search Area North Coordinate 

500 < (Ox < 2500 

m 

5 < m;c < 25 

Search Area East Coordinate 

500 < CO y < 2500 

m 

5 <cOy <25 

North Coordinate 

0 < X < 3000 

m 

0<x <30 

East Coordinate 

0 < y < 3000 

m 

0 <~y < 30 

Heading 

— OO < ifj < oo 

rad 

— OO < ij/ < OO 

Turn Rate 

-0.3 < r < 0.3 

rad/s 

-30<f < 30 

Rudder Angle Input 

-0.5 <u< 0.5 

rad 

-0.5 <u< 0.5 

Velocity 

V = 2.5 

m/s 

V = 2.5 

Nomoto Gain Constant 

K = 0.5 

1/s 

K = 50 

Nomoto Time Constant 

T = 5.0 

s 

f = 0.05 

Poisson Scan Rate 

A = 0.2 

1/s 

A = 20 

Attenuation Coefficient 

a = 0.052 

dB/m 

a = 5.2 


It is important to verify that these constraints are, in fact, satisfied in the continuous domain 
as well. Moreover, we require optimal trajectories that can be implemented on autonomous 
vehicles. Therefore, as a practical consideration, we adopt the definition of feasibility used 
by Humi: 

Showing the feasibility of the generated solution can be done by control trajec¬ 
tory interpolation and state propagation using a Runge-Kutta algorithm. If the 
initial conditions and system dynamics can be propagated using the optimal 
control solution and it matches the [solver’s] generated trajectories, then the 
control solution is deemed feasible. [74] 

In other words, solutions with discrete trajectories {x(k), u(k)] are feasible if a vehicle can 
execute a smooth control trajectory u{t), interpolated through the solution’s u{k) nodes, 
and produce a state trajectory x(t) sufficiently close to the solution’s x{k) nodes. Planning 
algorithms can verify feasibility automatically, but the word “close” must be quantified 
first. Possible metrics for the similarity between two curves include the Frechet [132] or 
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HausdorfF [133] distance measures. The trajectory-planning algorithm proposed in [74], 
for example, performs automatic feasibility checks using a norm based on summing the 
Euclidean distances between the solution nodes and points along the propagated trajectory, 
evaluated at the solution nodes. In the following chapters, we will consider a solution to 
be feasible when its state-propagated trajectory a) does not violate problem constraints, 
and b) matches the solution trajectory when overlaid on a graphical plot. These criteria 
will verify that solutions obtained from a numeric solver are feasible, and also ensure that 
only feasible guesses are used to initialize the optimization. When necessary, e.g., for 
the automated analysis of inverse problems conducted in Chapter 5, we employ a numeric 
feasibility criteria similar to [74]. 

3.7 Initial Guess 

Most numeric optimization routines are initialized with an initial guess. For an optimal 
control problem, the guess is a candidate solution, complete with state and control tra¬ 
jectories {x{k), u{k)] at discrete time nodes. The solver evaluates the objective function 
using these trajectories. From there, it iteratively generates new candidate solutions that 
decrease the objective value, finally stopping its search when it reaches a local minimum. 
A good initial guess can influence the optimization by focusing the solver’s effort in smaller 
regions of the search space. As a result, initial guesses can dramatically reduce solution 
times [74]. In some cases, e.g., when a problem has several local minima, the initial guess 
can determine whether a solver succeeds or fails at finding the correct solution. 

The initial guess must be a valid candidate solution to the optimization problem, which 
implies that the initial guess should: 

• have the same initial condition f (0) = xq as the problem of interest, 

• have the same time node discretization as the problem of interest, and 

• be feasible, i.e., obey state variable constraints and control limits. 

The first two requirements are easy to address while encoding the problem of interest into 
the GenOC framework. Satisfying the feasibility requirement depends on the sophistica¬ 
tion of the initial guess trajectory, which corresponds to the amount of prior information 
we wish to incorporate. Ideally, we would like to find an optimal solution without know¬ 
ing beforehand what a “good” initial guess looks like. A trivial guess which satisfies the 
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Figure 3.2. Initial guess trajectories for open loop rudder step. 

first two requirements, for example, would be a zero velocity trajectory that remains at xq 
for all time. Unfortunately, the constant velocity vehicle models defined in Section 3.1 do 
not permit acceleration, and the solver would be unable to find another solution trajectory. 
Similarly, a guess which specifies a trivial control trajectory u(0 = 0 for all time is infeasi¬ 
ble under our definition, because the vehicle would travel at constant velocity and heading 
until it left the domain of its North or East coordinate. 

A naive open-loop control trajectory is a good compromise between a trivial (no infor¬ 
mation provided) guess, and an expert (full information provided) guess. For example, 
a rudder angle step function, executed at the proper time, will cause a search vehicle to 
turn in a circle until the end of the simulation. This has the benefit of keeping the vehicle 
in the search area and ensures that state variable limits are not exceeded. In practice, we 
approximate a discrete step function with a smooth sigmoidal curve centered at the step 
time [96]. This simple control trajectory is then propagated through the motion model, us¬ 
ing a Runge-Kutta algorithm (e.g., the Matlab ode45 solver) to calculate the corresponding 
state variable trajectories, thereby guaranteeing feasibility (see Figure 3.2). 

If there is sufficient time to exhaustively search an area, an expert initial guess can be 
provided that completely covers the area with a deterministic search pattern. A num- 
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ber of algorithms for “coverage path planning” exist [63], based on sensor sweep width. 
Deterministic search patterns include spirals for searching circular areas [72]; and box- 
spirals, lawnmower, or zamboni patterns for searching rectangular areas [134]. While these 
strategies may waste effort when the search area contains sub-regions with near-zero tar¬ 
get probability distribution [135], it is usually possible to decompose the search area into 
smaller regions and avoid this situation. Moreover, for rectangular search areas, line sweeps 
conducted parallel to boundary edges are optimal for minimizing the number of turns re¬ 
quired [66]. This fact, plus the ease of implementing these patterns with vehicle autopilots, 
likely explain the widespread use of lawnmower patterns for underwater search operations. 

These rectangular coverage patterns require path-following controllers to execute them. 
Coverage path planners often take this for granted, assuming that the search area can be 
decomposed into smaller, “easy to cover” cells; the vehicle need only visit all such cells to 
achieve complete coverage [63]. Another approach is to extend the line sweep track length 
by a vehicle-specific distance, assuming all 180-degree turns occur outside the search area 
and the vehicle re-establishes straight-line motion before reentry on an adjacent track [136]. 
While these approaches determine the geometric length, spacing, and number of track lines 
for a given sweep width, they do not represent feasible trajectories, per se. 

Before we can specify this type of coverage pattern as an initial guess for an optimal control 
problem, we need to convert a waypoint-based specification into a feasible trajectory. Hau¬ 
gen suggests an approach for constructing a feasible lawnmower path which uses clothoids 
as transition curves between waypoint segments. The clothoids are scaled such that a ve¬ 
hicle following this trajectory obeys feasibility constraints on its angular velocity and ac¬ 
celeration [137]. Depending on the track spacing and vehicle turning radius, one of three 
different U-turn paths are constructed to connect adjacent line sweep tracks. For our search 
problems, we assume that the lawnmower track spacing permits a piecewise U-turn com¬ 
prised of two 90-degree clothoids and a straight line segment (Case A) [137]. This produces 
a feasible path in the horizontal plane, which we can convert to a control trajectory u{k) by 
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inverting the seareher model of Seetion 3.1 as follows: 


Ax(k) 

Ay(k) 

iff(k) 


r(k) = 


r(k) = 


At(k) 
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u(k) 


x(k + 1) - x(k) 
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atan2(Ay(k), Ax(k)) 

^ Ax^{k) + Ay^{k) 

V 

ik(k + 1) - ik(k) 
At(k) 

r{k + 1) - r(k) 

At{k) 

^{r{k) + Tr{k)). 


(3.27) 


Figure 3.3 and Figure 3.4 illustrate lawnmower and box-spiral initial guesses, respeetively, 
eonstrueted using this elothoid method. Note that the inverse kinematie equations of Equa¬ 
tion (3.27) differentiate state trajeetories using the forward Euler method, whieh requires 
small, equally-spaeed time steps to ensure aeeuraey and feasibility of the derivatives. Small 
step size translates into a large number of nodes, whieh ean drastieally inerease solver run 
time. The latter makes a eompelling argument against supplying a deterministie seareh 
pattern for the initial guess. 
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Figure 3.3. Initial guess trajectories for open loop lawnmower pattern. 
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Figure 3.4. Initial guess trajectories for open loop box-spiral pattern. 
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CHAPTER 4: 

Application: Time-Limited Optimal Search 


As we have seen in Chapter 3, eomputational optimal seareh ean be tailored to model 
and solve many different MCM problems of interest. Optimal solutions, obtained through 
simulation, provide performanee benehmarks that ean inform mission planning under real- 
world resouree limitations. These resourees inelude the number and type of autonomous 
vehieles at the MCM eommander’s disposal, as well as the sensors these platforms ean 
earry. Typieally, however, the most important resouree is time. As Washburn notes in [23]: 

The great questions in seareh all involve time. We ask, “How long will detee- 
tion take?” or “What is the probability of deteetion in a fixed time?” Deteetion 
is inevitable, given suffieient time. The objeet of seareh planning is to speed 
things up. [23] 

Indeed, while a number of planning algorithms have been developed to aehieve eomplete 
eoverage of a seareh area (see, e.g., [63], [138]), most do not explieitly eonsider the ram- 
ifieations of time. Instead, time is a byproduet of the seareh vehiele’s veloeity and spatial 
trajeetory. A eommon metrie is the area eoverage rate, eomputed by multiplying a sen¬ 
sor’s nominal sweep width by platform veloeity. One example deseribed in [139] derives a 
lower bound for the time required by an aerial vehiele to follow a flight plan that aehieves 
eomplete sensor eoverage. This type of bound ean be informative when there is suffieient 
time to exeeute a given motion plan, but provides no guidanee for adjusting the plan if the 
bound exeeeds the allowable mission duration. 

When time is limited and eomplete eoverage is impossible, deterministie seareh patterns 
(e.g., lawnmower or box-spiral trajeetories) are faeed with two ehoiees: 1) exeeute the 
original motion plan as long as possible to aehieve 100% sensor eoverage in a subset of 
the entire seareh area; or 2) adjust the traek spaeing to survey the entire seareh area, but 
with ineomplete eoverage. For a perfeet “eookie eutter” sensor and uniform target PDF, 
both ehoiees are equivalent and the probability of non-deteetion equals the fraetion of un- 
searehed area. Time-limited MCM operations ean only reduee this risk by leveraging prior 
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information about the target distribution. If it is known, for example, that mines have 
been deployed in an “evenly-spaced mine line,” [140] proposes a track-spacing method 
that yields a probability of missed mines below the un-searched area ratio. 

Mission time is a hard constraint in most MCM operations, motivating the fixed-time prob¬ 
lem formulation described in Chapter 3. We seek time-limited optimal search trajectories 
that minimize MCM risk for a given vehicle, sensor, and mission duration—whether or not 
prior information is available. In this chapter, we demonstrate the flexibility of the GenOC 
framework by solving MCM search problems for both cases. First, we implement an RID 
mission in which a vehicle with high-resolution sonar must revisit a previously detected 
mine target whose location is given by a joint, normalized, beta distribution. Next, we 
implement a wide area survey to detect and localize MLOs, using the conservative assump¬ 
tion of a uniform target PDF to represent no prior information. Such missions are typically 
conducted during the initial phase of an MCM operation. For this mission, we compare 
time-limited search performance, i.e., the MCM risk vs. allotted mission time, for optimal 
trajectories and well-known deterministic search patterns. Finally, we discuss the impact 
of time discretization on numerical solutions. 


4.1 Search with Prior Information—Mine Reacquisition 

During the initial phase of an MCM operation, wide area surveys are conducted to detect 
mine-like objects (MLOs) in the environment that pose a threat to naval forces. These 
surveys can produce datasets with dozens of potential target locations, so it is critical to 
distinguish actual mines from harmless clutter before launching time-intensive neutraliza¬ 
tion missions. Successful target identification requires high-resolution sensors not typically 
carried on the initial survey vehicles, and these sensors are more effective at close range. As 
a result, follow-on missions are conducted to revisit the MLOs with AUVs carrying imag¬ 
ing sonars or video cameras. This type of reacquire-identify (RID) mission incorporates 
prior information about MLO locations provided by the survey, but this data is uncertain; 
its accuracy depends upon the sensing and navigational performance of the survey vehicle 
itself. The searcher should expect to search for the target once it arrives in the vicinity of 
the surveyed location. We therefore cast the motion planning problem for an RID mission 
as an optimal search for a target whose probability density is more informative than the 
uniform density assumed for the initial MCM search. 
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In this section, we apply the GenOC framework to solve problems of this type. We assume 
that a prior survey has detected and localized a MLO with probability density described 
by a joint normalized beta distribution in two dimensions (see Section 3.3). For a search 
area of 20 DU x 20 DU in size, we select {a,P) parameter values of [8 DU, 16 DU] in 
the north direction and [16 DU, 8 DU] in the east direction. Substituting these values into 
Equation (3.6) produces the two dimensional PDF shown in Figure 4.1, pictured prior to 
commencement of the RID mission. Our objective value, the residual risk of non detection, 
is plotted on a color scale in which high probabilities are shown in dark red, and low 
probabilities are shown in blue. This jS distribution corresponds to a previously-detected 
target located at cOy] = [12 DU, 19 DU]. 

We wish to compute the time-limited optimal trajectory for a 40-minute RID mission by 
a REMUS 100 AUV equipped with high-resolution imaging sonar. For this problem, we 
assume that the search area shown in Figure 4.1 has a flat bottom and the AUV oper¬ 
ates with constant velocity V = 1.5 m/s at altitude h = 3 meters above the sea floor. 
The AUV is programmed to launch from a start location dX{x,y] = [1 DU, 7 DU] on 
an initial heading of 45 degrees, utilizing prior information about the expected target lo¬ 
cation [12 DU, 19 DU]. The search vehicle’s initial state vector is therefore given by 
T(0) = [1 DU, 7 DU, ;7r/4 rad, 0 rad/TU]^. A naive, yet feasible initial guess trajectory 
is provided to the solver using an open loop rudder step function to generate a wide right 
turn. 

The search performance of the 450 kHz and 900 kHz blazed array FLS models from Ta¬ 
ble 2.1 are compared in Figure 4.2 and Figure 4.3, respectively, for different numbers of 
discrete time nodes. To compare optimal trajectories generated with different discretization 
grids (in the time and/or spatial domains), we re-compute objective values using a common 
baseline of 500 time nodes and 25 x 25 grid of spatial nodes. 
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Figure 4.1. Normalized beta distribution for a previously detected target 
located at [12 DU, 19 DU]. 


We accomplish this by: 

1. Interpolating the solver’s control input trajectory onto a hne grid of 500 time nodes. 

2. Propagating this hnely-gridded control input through the vehicle’s EOM, using 
MATLAB’s ode45 solver to generate hnely gridded state trajectories. 

3. Re-calculating the objective function for the resulting 500-node, dynamically feasi¬ 
ble state trajectories. 
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Figure 4.2. Optimal RID trajectories for a REMUS AUV with P450 FLS. 
Left: Pnd = 0.012 (20 time nodes). Right: Pnd = 0.001 (35 time nodes). 



Figure 4.3. Optimal RID trajectories for a REMUS AUV with P900 FLS. 
Left: Pnd = 0.333 (30 time nodes). Right: Pnd = 0.305 (40 time nodes). 


The 20-node trajeetory in Figure 4.2 (at left) makes two parallel passes over the seareh area 
and then loiters over the expeeted target loeation for the remainder of the mission. The 35- 
node trajeetory shown at right is more dynamie, approaehing the seareh area from several 
different headings and almost encireling the expeeted target location. The 450 kHz blazed 
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array FLS has a nominal operating range of 200 meters, and an AUV using this sensor for 
RID missions can reduce the residual risk of non-detection to roughly 1% or less. 


If a given mission requires a higher-resolution sonar to achieve positive target identification, 
the same AUV can deploy a 900 kHz blazed array FLS instead. Example trajectories are 
provided in Figure 4.3. The 30-node trajectory shown at left executes three symmetric loops 
over the expected target location, while the 40-node trajectory shown at right sweeps the 
target area with parallel tracks oriented on two distinct headings. Both solutions achieve a 
probability of non-detection around 30%, indicating that the 900 kHz sonar is hampered (in 
coverage) by its 100-meter nominal range. Nevertheless, the GenOC framework produces 
trajectories that revisit the target area repeatedly until the mission time expires. These 
solutions have similarities with traditional RID search patterns, which implement partial 
lawnmower swaths, aligned on different headings, to cover a target from multiple aspect 
angles for improved classification performance [30]. Results from these simulations are 
summarized in Table 4.1. 

Table 4.1. Optimal time-limited RID trajectories for a REMUS 100 AUV 

with FLS and Tf = 2400 s. 


Imaging 

Time 

Pnd 

Ligure 

Sonar 

Nodes 

500 X 25 X 25 

Reference 

P450 FLS 

20 

0.012 

Ligure 4.2 (left) 

35 

0.001 

Ligure 4.2 (right) 

P900 FLS 

30 

0.333 

Ligure 4.3 (left) 

40 

0.305 

Ligure 4.3 (right) 


4.2 Search with No Prior Information—Mine Survey 

In this section, we consider the common MCM problem of planning survey missions to de¬ 
tect MLOs in the absence of prior information about the target distribution. Typically, this 
involves a labor-intensive process to divide the search area into separate homogeneous re¬ 
gions, and tasking individual MCM assets to cover each region with a deterministic search 
pattern based on nominal area coverage rates. Various TDAs have been developed to au¬ 
tomate aspects of this process, and these information products can be “used by the force 
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commanders to optimize the employment of naval assets in any partieular taetieal environ¬ 
ment at sea” [56]. 

Many planning tools, e.g., PATHA, ean ineorporate sonar performanee models and time- 
based eonstraints to help mission planners determine the number of assets needed for a 
given mission [59]. However, these systems do not explieitly eonsider search vehiele dy- 
namies and their attendant impaet on detention performanee. The GenOC framework takes 
this into aeeount, providing a unique eapability for planning MCM survey operations. It not 
only produees time-limited optimal seareh trajeetories that minimize risk for a given vehi¬ 
ele and sensor eonfiguration, but it ean also serve as a pre-mission analysis tool to eompare 
the expeeted performanee of different sonar designs or autonomous vehiele teams. The lat¬ 
ter refers to the solution of so-ealled “inverse problems,” whieh is diseussed in Chapter 5. 
The rest of this ehapter is foeused on the former, and applies the GenOC framework to plan 
a time-limited MCM survey for the following benehmark problem. 

We wish to plan 30-minute MCM survey that aehieves 90% probability of deteeting a 
bottom mine hidden anywhere in the 20 DU x 20 DU seareh area shown in Figure 4.4, 
with uniform probability distribution. This eorresponds to a risk threshold ^ 10%. 
For this problem, we assume that the seareh area has a flat bottom, water depth is 20 meters 
(e.g., 0.2 DU), and two SeaFox US Vs equipped with 200 kHz FLS are available for this 
operation. We first eonsider whether a single vehiele ean meet the desired risk threshold, 
and solve an optimal search problem for a SeaFox USV launehed from the initial state 
veetor f(0) = [1 DU, 7 DU, 0 rad, 0 rad/TU]^ and programmed to operate at eonstant 
veloeity V = 2.5 m/s = 2.5 VU. A naive, yet feasible initial guess trajeetory is provided to 
the solver using an open loop rudder step funetion to generate a right turn in the eenter of 
the seareh area. 
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East [100m] 


Figure 4.4. Search area with uniform probability distribution, representing 
no prior target location data. 


The optimal time-limited seareh trajeetory is shown in Figure 4.5, whieh aehieves an 
objeetive value of Pnd = 0.362, eomputed on our 500 x 25 x 25 diseretization base¬ 
line. This trajeetory resembles a box-spiral seareh pattern, although the limited mission 
duration does not permit full eoverage. Nevertheless, this trajeetory represents the best 
seareh performanee that ean be aehieved by a single USV launehed from the given ini¬ 
tial eondition. As sueh, it represents a loeal minimum, sinee different initial eonditions 
may yield lower Pnd results. Monte Carlo simulation ean be employed to determine the 
most favorable initial eondition(s) for exploring a given seareh area. Due to the sym¬ 
metry of this problem, however, sueh solutions are not unique. An initial eondition of 
x(0) = [23 DU, 1 DU, 7r/2 rad, 0 rad/TU]^ will yield the same result. 
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Figure 4.5. Optimal survey trajectory to detect MLOs with a SeaFox USV 
and 200 kFIz FLS in Tf = 30 minutes: P^d = 0.362 (30 time nodes). 

Since the MCM survey failed to meet the desired risk threshold, the foree eommander must 
either inerease the mission duration, or deploy additional search assets. Adding a seeond, 
identical searcher launehed from the initial eondition x(0) = [1 DU, 9 DU, 0 rad, 0 rad/TU]^ 
produees the optimal trajeetories shown in Figure 4.6, redueing the risk to = 0.022 in 
the same 30-minute mission. This result suggests that two vehieles ean meet the desired 
risk threshold of P^d =0.1 with a shorter mission duration. 
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Figure 4.6. Optimal survey trajectories to detect MLOs with two SeaFox 
USVs and 200 khlz FLS in Tf = 30 minutes: = 0.022 (30 time nodes). 


Additional simulations were conducted to determine whether two vehicles could achieve 
the survey objective in less time. Figure 4.7 shows the optimal trajectories computed for 
a 20-minute (at left) and 25-minute (at right) mission. The 20-minute mission has time- 
limited performance = 0.224 and fails to meet our objective, while the 25-minute 
mission achieves Pj^d = 0.087, equivalent to a detection probability Pd > 91%. Results 
from these simulations are summarized in Table 4.2. 
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Figure 4.7. Optimal mine detection survey trajectories for two SeaFox USVs 
with 200 kHz FLS. Left: Pnd = 0.224 (20 time nodes) in 20 minutes. 
Right: Pi^D = 0.087 (30 time nodes) in 25 minutes. 


Table 4.2. Optimal time-limited trajectories for mine detection using SeaFox 
USVs and 200 kHz FLS. 


Number 

Mission 

Time 

Pnd 

Figure 

of USVs 

Duration 

Nodes 

500 X 25 X 25 

Reference 

1 

1800 s 

30 

0.362 

Figure 4.5 


1200 s 

20 

0.224 

Figure 4.7 (left) 

2 

1500 s 

30 

0.087 

Figure 4.7 (right) 


1800 s 

30 

0.022 

Figure 4.6 


4.3 Search Performance vs. Mission Duration 

There is an inherent time-dependency of the exponential detection model incorporated into 
the objective function of Equation (3.11). This model produces diminishing returns on 
search effort applied to previously visited regions of the operating area. Optimal search can 
leverage this property to produce motion plans which accomplish both exploration, when 
we wish to acquire information about the environment (Section 4.2); and exploitation of 
all relevant prior information (Section 4.1). We have already seen how mission duration 
impacts the optimal vehicle trajectories and achievable search performance for a given 
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mission. While it is intuitively obvious that searching for longer periods of time can lower 
MCM risk, the ability to rapidly solve optimal search problems allows MCM commanders 
to quantitatively address questions such as: 

• What is the residual risk after searching for a fixed duration with a given vehicle and 
sensor combination? 

• How long will it take a search vehicle to reach a desired risk threshold with a given 
sensor payload? 

• How much time savings can be achieved by employing multiple search assets? 

Using the GenOC framework, we can conduct several simulated experiments to solve a 
given optimal search problem for different values of our fixed mission duration. Subsequent 
Monte Carlo analysis can identify trends in the solution results to help characterize the 
optimal performance of a given system configuration as a function of time. Note that our 
current computational framework was written to solve GenOC problems with fixed final 
time; at present, we rely on Monte Carlo simulations to answer questions regarding the 
minimum time to reach a given risk threshold, for example. Future work will investigate the 
use of optimal control software packages l ik e DIDO [141], which can address minimum¬ 
time problems directly, for solving this class of GenOC problems. 

Meanwhile, we must address two minor complications arising from this approach. First, 
the numeric solution to a given optimal control problem is repeatable over multiple runs; we 
must therefore inject random variation into our simulations before we can conduct a mean¬ 
ingful statistical analysis. This is achieved by varying the initial condition; the position and 
heading angle of each search vehicle at t = 0 is randomized prior to each simulation, influ¬ 
encing the initial guess trajectories as well. Second, we must impose a feasibility check on 
the solver’s output so that infeasible solutions can be excluded from the analysis. 

It is instructive to compare the search performance vs. mission duration of our optimal 
search trajectories against well-known deterministic search patterns. For the single-USV 
survey mission described in Section 4.2, the following sections describe the feasible lawn- 
mower and box-spiral trajectories which provide benchmarks for comparison. 
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4.3.1 Manually-Specified Lawnmower Pattern 

Computing the shortest path for a lawnmower coverage pattern to completely cover a 
polygonal area has been found to be NP-hard; as a result, approximate algorithms are pro¬ 
posed to plan efficient lawnmower trajectories in [142]. Additional examples are proposed 
in [66], which minimizes the number of turns along the path, and [143] which suggests that 
these patterns are time-optimal for a robotic lawnmower. 

We wish to generate a benchmark lawnmower trajectory that completely covers the 
20 DU X 20 DU search area shown in Figure 4.4. We do so in the manner that an MCM 
operator programs a waypoint-based mission for a given sensor sweep width. The 200 kHz 
FLS used in our example has a nominal range of 400 meters, so we select waypoints that 
place north/south-aligned track lines with 400 meters track spacing, offset by 200 meters 
from the search area’s east/west boundaries. Whereas a sidescan sonar survey would per¬ 
form the necessary U-turns outside the search area, we assume that the FLS does not incur 
a performance penalty for turning motions (i.e., = 1). Therefore, the long track lines are 

connected by short legs offset 200 meters from the search area’s north/south boundaries to 
avoid wasted effort. This waypoint pattern is denoted by green circles in Figure 4.8. 

Similarly, we specify an initial condition of v(0) = [1 DU, 7 DU, 0 rad, 0 rad/TU]^, so the 
USV begins its mission pre-aligned with the first track line, thereby minimizing unneces¬ 
sary path length. Finally, we ensure feasibility of this lawnmower pattern by constructing 
each 90-degree turn with clothoid curves (see Section 3.7) and propagating a control input 
through the EOM to produce the magenta state trajectory shown in Figure 4.8. This trajec¬ 
tory was used to recalculate Pnd on our 500 x 25 x 25 discretization baseline for different 
values ofTf. These results are plotted as blue crosses in Figure 4.9, where each data point 
represents an entire mission that completes as much of the lawnmower pattern as possible 
in the time allotted. 
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Time = 20 [100s] 




Figure 4.8. Partially completed lawnmower survey for = 2000 s. 


Figure 4.9 reveals a linear relationship. We note that P^d falls at a constant rate for mission 
durations up to roughly 3500 seconds, then slows down just before vanishing altogether 
for missions around 4000 seconds long. The best linear fit to this data is obtained by 
eliminating data points when 7/ > 3500 seconds, producing the formula 

Pnd = -0.000253{Tf) + 0.950. (4.1) 

Equation (4.1) confirms what we would expect: Pnd ~ 1 for Tf = 0 seconds, a trivial 
mission with no search effort; and Pnd = 0 for Ty > 3750 seconds, which is sufficient time 
to completely cover the search area. 
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MCM Risk vs. for Deterministic Search Patterns 


(1 USV searching 4 km^ area) 
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Figure 4.9. Linear fit of lawnmower and box-spiral search performance vs. 
mission duration. 


4.3.2 Manually-Specified Box-Spiral Pattern 

Spirals are another popular deterministic search pattern, and box-spirals are very similar 
to the lawnmower trajectories discussed previously. While both patterns achieve complete 
coverage, [143] suggests box-spirals as a minimal-energy alternative to time-optimal lawn- 
mower patterns, since box-spirals require less turning effort. Using the same assumptions, 
waypoint spacing, and initial condition as Section 4.3.1, we generate a benchmark box- 
spiral trajectory that completely covers the 20 DU x 20 DU search area of Figure 4.4. 
The feasible state trajectory for this pattern is the magenta line in Figure 4.10, while the 
waypoint pattern is denoted by green circles. 
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Figure 4.10. Partially completed box-spiral survey pattern for Tf = 2800 s. 


Pnd values, recalculated on our 500 x 25 x 25 discretization baseline, are plotted as red 
boxes in Figure 4.9 for different values ofTf. As before, each data point represents an entire 
mission that completes as much of the box-spiral pattern as possible in the time allotted. 
This data reveals a linear relationship nearly identical to the lawnmower pattern. The best 
linear fit yields the formula: 


Pnd = -0.000254(r/) + 0.947. (4.2) 

The mean of the lawnmower and box-spiral linear fits is depicted by a black line with 
constant slope in Figure 4.9. A constant slope is expected, since lawnmower and box-spiral 
trajectories implement “exhaustive search.” Recall that exhaustive search with a definite 
range sensor yields a detection probability Pd equal to the coverage ratio, the fraction of 
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search area covered by the sensor. Therefore, a “cookie cutter” sensor with sweep width 
W, mounted on a vehicle moving at constant velocity V, in a search area A, produces a 
probability of detection Pd = WVt/A that is linear with time [51]. We have simply plotted 
the complement, Pnd = 1 - Pd, in Figure 4.9. 

4.3.3 Optimal Solution from Solver 

Monte Carlo simulations were conducted for the mine survey problem of Section 4.2 to 
analyze the optimal search performance for a range of mission durations between fifteen 
minutes and one hour, spaced at two-minute intervals. Ten simulations were conducted 
for each value of Tf, with initial states drawn from a uniform 'U (lower, upper) or normal 
AfCmean, std. dev.) probability distribution as follows: 

x{0) - 'I/(0,10)/DC, 

i/(0) ~ nU{500,2500)1DU, and (4.3) 

i/^(0) ~ Af(0,7r/12). 

Initial turn rate r(0) = 0. The initial guess is computed from an open loop rudder step 
input that commands a right turn when y (0) is in the western half of the search area, and 
commands a left turn when y (0) is in the eastern half of the search area. 

The mean objective values of Pnd from the ten simulations conducted for each value oiTf 
are plotted as the blue line in Figure 4.11, while a quadratic curve fit to this data is shown 
as the dashed magenta line. The mean of the two deterministic search patterns derived 
previously is shown as a black line for comparison purposes. The plots intersect at roughly 
Tf = 3280 seconds, at an objective value of Pnd = 10%. This plot clearly indicates that 
an optimal search strategy outperforms deterministic, exhaustive search patterns for time- 
limited missions less than about 55 minutes in duration. However, it also suggests that 
operators would be better off selecting a deterministic search pattern if there is sufficient 
time to execute it to completion. 
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MCM Risk vs. for Deterministic Search Patterns 


(1 USV searching 4 km^ area) 



Figure 4.11. Search performance comparison between optimal trajectories 
and exhaustive search patterns. 


The cross-over point may indicate a local minimum that is difficult to overcome with tra¬ 
jectories based on a low number of discrete time nodes. Simulations with excess mission 
duration are capable of bringing Pj^d to zero (see RID trajectories in Figure 4.2), albeit 
with some wasted search effort. In addition, previous simulations conducted using an Euler 
discretization with hundreds or thousands of time nodes have generated trajectories that 
more closely resemble lawnmower patterns, but this kind of overkill is computationally 
prohibitive and not useful for optimal trajectory planning. Futher work is required to inves¬ 
tigate optimal search performance when mission duration is roughly equivalent to the time 
required for exhaustive search. 


98 





































4.4 Search Performance vs. Time Discretization 

One of the most important aspeets of motion planning in a eomputational optimal eontrol 
framework is the ehoiee of diseretization seheme, as this direetly impaets aeeuraey and 
eomputational run time. In general, numerie trajeetory approximations eonverge to their 
eontinuous eounterparts as the number of eomputational nodes inerease [144], [145]. While 
inereasing the number of nodes ean improve solution aeeuraey, designers must balanee 
this aeeuraey against the eomputational demands required by high-node diseretizations. 
Moreover, higher-node eontrol trajeetory solutions may be infeasible for implementation 
on an aetual vehiele system. A detailed theoretieal diseussion on this topie is beyond the 
seope of this dissertation, but an exeellent overview on pseudospeetral optimal eontrol 
theory is provided in [4], with eonvergenee and eonsisteney proofs given in [77]. 

Hurni reeommends using “the lowest possible number of nodes for feasible and safe tra- 
jeetories,” and proposes a novel eriteria for seleeting the number of nodes based on the 
distanee a ground vehiele must travel, and the size of the obstaeles it must avoid along 
the way [74]. We have assumed an obstaele-free environment for MCM seareh planning. 
Moreover, we do not require real-time algorithms for dynamie re-planning. Instead, we 
generate an optimal seareh strategy for the entire MCM mission, subjeet to any prior infor¬ 
mation we possess. The laek of a real-time eonstraint grants us the luxury of eomputing 
multiple solutions with inereasingly fine diseretization sehemes in our seareh for an opti¬ 
mal, feasible seareh trajeetory. Solutions whose feasibility eannot be verified by eontrol 
trajeetory propagation are rejeeted. 

We will demonstrate this eoneept using the time-limited survey mission deseribed in See- 
tion 4.2. Figure 4.12 plots the objeetive values for ten different optimal solutions, eaeh 
eomputed using a different time diseretization with Nt time nodes. We observe that as 
the number of time nodes inerease, the numerie solution’s objeetive value eonverges to 
Pnd ~ 0.320. Note, however, that the objeetive values plotted in Figure 4.12 are the raw 
solver outputs, ealeulated direetly from the Nt solution trajeetory nodes. We denote sueh 
objeetive values by Jout- 

To support performanee eomparisons between different numerie solutions, the solution tra- 
jeetories must first be transferred onto a eommon diseretization seheme. This is aehieved by 
interpolating eontrol trajeetories, propagating state variables through the system’s ODEs, 
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and recalculating objective values on our 500 x 25 x 25 discretization baseline. Objective 
values corresponding to these ODE-propagated state trajectories are designated by Jode- 
Finally, the propagated trajectories are compared against the solution nodes to assess each 
solution’s feasibility. In this example, all of the numeric solutions with greater than 50 
nodes are deemed infeasible. Feasible search trajectories are illustrated in Figure 4.13. 

Time-limited Performance vs. Time Nodes 


Search with 1 USV for Tf = 18 [TU] 



Figure 4.12. Solver-provided objective values vs. the number of discrete 
time nodes for a single-vehicle, 30-minute survey mission. 
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Optimal Trajectories vs. Time Nodes for Tf = 18 [TU] 



Figure 4.13. Single-vehicle, 30-minute survey trajectories computed for dif¬ 
ferent time discretizations. 


Note that higher-node solutions ineorporate periodie turning motions in their trajectories. 
This has the benefit of aiming the vehicle’s FLS to cover a larger portion of the search area, 
reducing the accumulated probability of non-detection. Recall that we have not penalized 
the detection rate of this FLS for turning motion as we would for a SSS, i.e., the shaping 
function Fr{x) = 1 for this problem. This increased complexity yields diminishing returns, 
however; the optimal trajectories computed using 30 or more time nodes are remarkably 
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similar, and all of them achieve a Pj^id within 3% of the 50-node best performer. 

Figure 4.14 illustrates the run times required to compute each solution in Figure 4.12 and 
Figure 4.13 using a 2.30 GHz Xeon CPU (complete processor specifications are listed 
in [146]). We note that there is a large increase in run time required to compute solutions 
with more than 40 time nodes, and the 50-node solution takes nearly twice as long as the 
30-node solution. Moreover, the 40-node solution takes nearly 5 seconds longer to compute 
than the 30-node solution, but only decreases by 0.003. 


CPU Run Time vs. Time Nodes 
Search with 1 USV for Tf = 18 [TU] 



Figure 4.14. Run times vs. the number of discrete time nodes required to 
compute trajectories for a single-vehicle, 30-minute survey mission. 


For this reason, 30-node solutions were selected as an acceptable compromise for the plots 
presented in Section 4.2. Table 4.3 summarizes the search performance for several one- and 
two-vehicle mine detection surveys, computed for different mission durations with different 
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numbers of discrete time nodes. The corresponding two-vehicle trajectories are shown in 
Figure 4.15 through Figure 4.17. 


Table 4.3. Optimal time-limited search performance vs. the number of 
discrete time nodes for mine detection survey missions. 


No. of 

USVs 

Mission 

Duration 

Time 

Nodes 

Jout — Pnd 

NX 25x25 

Jode — P ND 

500 X 25 X 25 

CPU 

Run Time 

Figure 

Reference 

1 

1800 s 

10 

0.676 

0.431 

14.1 s 

Figure 4.13 

20 

0.443 

0.353 

28.0 s 

30 

0.358 

0.337 

31.4 s 

40 

0.342 

0.334 

36.0 s 

50 

0.328 

0.328 

63.1 s 

2 

1200 s 

10 

0.453 

0.303 

47.8 s 

Figure 4.15 

20 

0.218 

0.188 

104.7 s 

1500 s 

10 

0.428 

0.251 

47.1 s 

Figure 4.16 

20 

0.160 

0.133 

77.5 s 

30 

0.076 

0.067 

124.5 s 

40 

0.058 

0.056 

141.5 s 



10 

0.459 

0.273 

46.3 s 



1800 s 

20 

0.117 

0.065 

83.7 s 

Figure 4.17 



30 

0.021 

0.011 

119.5 s 
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North [DU] 


Optimal Trajectories vs. Time Nodes for Tf = 12 [TU] 



Figure 4.15. Two-vehicle, 20-minute survey trajectories computed for differ¬ 
ent time discretizations. 
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North [DU] 


Optimal Trajectories vs. Time Nodes for Tf = 15 [TU] 



Figure 4.16. Two-vehicle, 25-minute survey trajectories computed for differ¬ 
ent time discretizations. 
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North [DU] 


Optimal Trajectories vs. Time Nodes for Tf = 18 [TU] 



Figure 4.17. Two-vehicle, 30-minute survey trajectories computed for differ¬ 
ent time discretizations. 
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CHAPTER 5: 

Application: Inverse Problems 


The ability to rapidly solve optimal search problems provides a new tool for investigating 
so-called “inverse” problems related to optimal vehicle and sensor configurations for mine 
countermeasures (MCM). In Chapter 4, we used the GenOC framework to answer “direct” 
questions related to sensor-based motion planning for MCM. Typical direct applications 
involve calculating a feasible trajectory that minimizes the risk of non-detection at the 
conclusion a wide area survey using a specific sonar; or determining the time required 
for an autonomous vehicle team to achieve a given performance threshold. Conversely, 
inverse problems provide engineering insights about the search assets themselves. Such 
questions include: What is the most effective sonar mounting angle for an optimal search 
conducted by USVs? Which sonar design parameters have the biggest impact on search 
performance? How do multiple low-cost systems perform in a given mission, compared 
against the performance of a single expensive asset? 

Over the years, acoustic modeling and simulation have played an increasingly valuable role 
in the design and evaluation of sonar systems for naval operations. Etter observes: 

Modeling has become the chief mechanism by which researchers and analysts 
can simulate sonar performance under laboratory conditions. Modeling pro¬ 
vides an efficient means to parametrically investigate the performance of hy¬ 
pothetical sonar designs ... and to estimate the performance of existing sonars 
in different ocean areas and seasons. [56] 

Indeed, the Generic Sonar Model, a predecessor to the USN’s current CASS system, was 
originally “designed to provide sonar system developers with a comprehensive modeling 
capability for evaluating the performance of sonar systems and ... permit cost/accuracy 
trade-offs for specific applications” [147]. 

This chapter describes a novel approach to solving these kinds of inverse problems, one 
which recognizes the fact that sonar effectiveness often depends on the vehicles that deploy 
them. The physics-based sonar models derived in Chapter 2 provide a link between search 
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Figure 5.1. Analysis of multi-vehicle search performance with different 
sonars. Each data point represents several optimal vehicle trajectories. 


performance and sonar design parameters, while computational optimal search facilitates 
the rapid solution of multiple problems required for Monte Carlo analysis. In this way, 
sonar designers and operational planners can numerically determine the sensitivity of a 
given optimal search scenario to individual vehicle and/or sonar design parameters. Of 
course, a beneficial by-product of this analysis is a set of optimal vehicle trajectories linked 
to a given scenario. This is depicted in Figure 5.1, which emphasizes that each data point 
in a Monte Carlo analysis is comprised of several optimal trajectories that produced that 
result. 
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The Hamming supercomputer, a “hybrid cluster” of computing cores available via the NFS 
High Performance Computing (HPC) Center [148] was used extensively during the course 
of this dissertation. Although our computational framework utilizes a sequential optimiza¬ 
tion algorithm that doesn’t lend itself to massive parallelization, several simulations can be 
launched via the network to run on multiple computing cores simultaneously. This capabil¬ 
ity was a key enabler for conducting numerous trade studies at once, greatly accelerating 
the analysis process. Copies of the batch files used to generate the Monte Carlo simulation 
results for this chapter are provided in Appendix B. 

5.1 Sensitivity Analysis 

Sensitivity is used to describe the manner in which a desired output variable changes in re¬ 
lation to a change in some other system parameter. Sensitivity analysis provides a measure 
of robustness in control system design. For our optimal search problem, the desired output 
is our objective value, the residual MCM risk upon completion of a search operation. To 
address so-called inverse problems, it is helpful to calculate the sensitivity of MCM risk as 
a function of sonar/vehicle design parameters or operating characteristics. While our objec¬ 
tive function can be expressed analytically, deriving an analytic expression for sensitivity 
is difficult; it is easier to compute this metric numerically, especially after multiple Monte 
Carlo simulations have produced ample numerical data to work with. This is the approach 
taken in this chapter. We note, however, that optimal control offers an elegant method for 
computing sensitivity via parameter value-function analysis. The Lagrange multiplier for 
a parameter of interest in an optimal control problem serves as a “sensitivity multiplier.” 
Therefore, it accurately describes the relative impact of this parameter on the objective 
function [149]. The DIDO software package computes these costates during the solution 
of an optimal control problem and provides them as outputs, greatly facilitating this type 
of sensitivity analysis. 

5.2 Objective Function Value vs. Time Discretization 

As we observed in Section 4.4, the final objective value for optimal solutions computed 
for an increasing number of time nodes eventually converges, but higher-node solutions are 
more likely to fail our feasibility criteria. Therefore, it is important to determine the number 
of time nodes which provides an accurate value for PNoiTf) yet still produces feasible 
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trajectories. We accomplish this through a Monte Carlo analysis of several simulations 
conducted for an array of different time nodes. Once the acceptable number of time nodes 
has been identified for a given problem, subsequent simulations used to investigate other 
inverse problems can utilize the same discretization scheme. 

The parameter values used for this analysis are listed in Table 5.1, where the bold array no¬ 
tation indicates the range of the independent variable under investigation, i.e., simulations 
are performed for Nt = 20,25,30,... ,60. A separate analysis was conducted for each of 

Table 5.1. Simulation parameters for Nt analysis (free parameters in bold). 


Symbol 

Definition 

Value 

Tf 

Mission duration [min] 

30 

Nsims 

Number of simulation runs 

10 

Nt, 

Number of vehicles 

1 

Nt 

Number of time nodes 

[ 20 : 5 : 60 ] 

Nt, 

Number of parameter nodes 

25x25 

DU 

Canonical distance [m] 

100 

TU 

Canonical time [s] 

100 

UxV 

Vehicle dynamics file 

user-selected 

FLS 

Sonar parameters file 

user-selected 


the three FLS models derived in Chapter 2. We adopt feasibility criteria similar to the one 
proposed in [74] to reject infeasible solutions: 


r = 




i=l 


(5.1) 


where [xt, i/,] is the location of solution node i, and is the location along the prop¬ 

agated trajectory, interpolated at time t,-. We declare feasibility if < Tmax, where 
Tmax 6 [2,3] provides good rejection criteria for this problem. For each of the Nt val¬ 
ues listed in Table 5.1, Figure 5.2 plots the mean and standard deviation (depicted by solid 
lines and error bars, respectively) for the objective values from ten simulations. Figure 5.3 
shows the fraction of these search trajectories which meet our feasibility criteria. Search 
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performance for the solver output Jout and the ODE-propagated solution Jode (recalculated 
on our 500 x 25 x 25 discretization baseline) are both provided for comparison purposes. 


MCM Risk vs. Number of Time Nodes 
(1 USV searching 4 km^ area) 



Figure 5.2. Average single-vehicle search performance vs. number of nodes 

(N.y 


Fraction of Feasible Trajectories from 10 Time Node Simulations 



20 25 30 35 40 45 50 55 60 

Number of Time Nodes 


Figure 5.3. Fraction of single-vehicle Nt simulations with feasible trajectories. 
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These plots illustrate the trade-olF between accuracy and feasibility of higher-node solu¬ 
tions. In this problem, all of the 30-node solutions produced feasible search trajectories for 
all three sonar models, but their objective values have not yet converged to the objective 
values of their ODE-propagated trajectories. Alternatively, while the 55- and 60-node solu¬ 
tions yield objective values similar to their ODE-propagated trajectories, these higher-node 
solutions produce fewer feasible trajectories. On the basis of these results, therefore, we 
have elected to use 50 time nodes for the remaining Monte Carlo analyses presented in this 
chapter. 

5.3 Single-Vehicle Search Performance vs. Sonar Design 
Criteria 

We now consider several inverse problems which investigate how individual sonar design 
parameters influence optimal search performance for a given MCM search asset. This 
type of analysis can help sonar developers identify promising equipment modifications that 
could yield large performance benefits. It can also inform system operators about which 
configurable settings are most effective in a given scenario. Recall from the sonar design 
discussion in Chapter 2 that several different parameters influence a sonar’s detection rate. 
The operating frequency directly impacts the attenuation coefficient, which determines the 
propagation losses due to absorption. Erequency also plays a role in computing a sonar 
array’s directivity index and the ambient noise level the sonar must contend with. All of 
these influences enter the signal excess equations and impact a sonar’s instantaneous detec¬ 
tion probability. Other parameters, such as the Poisson Scan rate A contribute directly to 
the sonar’s instantaneous detection rate, the main driver of our objective function. Einally, 
geometric dependencies based on EOV or sonar mounting angle can make all the difference 
between an effective search operation, and one in which the residual risk of non detection 
is too great. 

We propose using the GenOC framework to asses the impact of sonar design parameters on 
search performance. The general approach used for the Monte Carlo analyses in each of the 
following sections is to hold all parameters constant except the parameter of interest, and 
conduct multiple optimal search simulations for each set of parameter values. We utilize 
50 time nodes for all simulations, based on the analysis in Section 5.2. The proposed 
approach requires certain choices to be made regarding which values to use for parameters 
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being held constant in a given simulation. In most cases (e.g., for FOM, Poisson Scan 
rate A, and horizontal FOV), we choose the median parameter value from among the three 
FLS systems. Moreover, we ensure a fair comparison among the different sonar designs by 
setting the horizontal field of view to 90 degrees for all systems. Lastly, for analyses not 
directly concerned with a given sonar parameter (e.g., optimal number of search assets), 
we utilize the actual sonar models listed in Table 2.1. 

An exception to this policy is made for the case of a sonar’s vertical mounting angle Vde, 
since initial simulation results that held Vde to a median value for all three sonars un¬ 
fairly penalized detection performance of the 200 kHz and 900 kHz sonars. In fact, vertical 
mounting angle Vde plays a critical role in a sonar’s effectiveness, particularly when search¬ 
ing for bottom mines with a USV. It is worth identifying other parameters which might be 
tightly-coupled through either problem geometry or sonar frequency. The FOM, for exam¬ 
ple, reflects a sonar design’s achievable range based on positive signal excess. We expect 
that this parameter would couple with the sonar’s vertical mounting angle, since these two 
parameters define a trigonometric relationship between the sonar and a bottom target. 

Fortunately, we have determined that using the median value of 66 dB for the FOM in 
certain simulations had only a minor impact on the problem’s geometry. Specifically, the 
noise-limited range for a 200 kHz FLS with nominal FOM = 72 dB decreases by about 
81 meters for a FOM of 66 dB. This only changes the optimal Vde from -5.4 to -6.1 degrees, 
a difference of only -0.7 degrees. Similarly, increasing the nominal FOM for the 900 kHz 
FLS from 64 dB to 66 dB increases noise-limited detection range by about 5 meters. This 
changes the optimal Vde for this sensor from -13.1 degrees to -12.4 degrees, a difference 
of only -1-0.7 degrees. Therefore, we justify the decision to analyze vertical mounting angle 
first, and use the optimal Vde values determined from these simulations in all subsequent 
analyses. 

5.3.1 Vertical Mounting Angle 

Vertical mounting angle is an important consideration for detecting bottom mines, as this 
parameter determines the ability of a sonar’s beams to cover the sea floor from a vehicle 
platform’s operating altitude. While sophisticated sonar like the ATLAS can electronically 
steer their beams in the vertical plane, lower-cost systems typically transmit at a fixed 
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angle. Therefore, these sonars are usually hard-mounted onto a vehiele at a vertieal angle 
that optimizes the sonar imagery eollected from the vehiele’s customary operating altitude. 
A custom BlueView P450 FLS system designed for use on the NPS REMUS 100 AUV, 
for example, was designed with multiple blazed arrays mounted at a permanent tilt angle 
of Vde = -10 degrees (see Figure 5.4). 




Figure 5.4. Custom BlueView P450 blazed array FLS mounted at -10 de¬ 
grees. 


The parameter values used for this analysis are listed in Table 5.2, where the bold array 
notation indicates the range of the independent variable under investigation, i.e., simula¬ 
tions are performed for Vde = -25,-20,-15,... ,-5 degrees. For each of the Vde values 
listed in Table 5.2, Figure 5.5 plots the mean and standard deviation (depicted by solid 
lines and error bars, respectively) for the objective values from ten simulations. Figure 5.6 
shows the fraction of these search trajectories which meet our feasibility criteria. Search 
performance for the solver output Jout and the ODE-propagated solution Jgde (recalculated 
on our 500 x 25 x 25 discretization baseline) are both provided for comparison purposes. 

Not surprisingly, longer-range sonars (i.e., with higher FOM) perform better at small 
mounting angles, while shorter-range, high-resolution sonars require steeper angles to ef¬ 
fectively reach the bottom. From this analysis, we can determine that the optimal mounting 
angles for detecting bottom mines in 20 meters of water with 200 kHz, 450 kHz, or 900 kHz 
FES from on a USV are -6 degrees, -7 degrees, and -11 degrees, respectively. Even so, the 
900 kHz FES is poorly suited for detecting bottom mines from a USV in this benchmark 
problem’s 20-meter water depth. 
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Table 5.2. Simulation parameters for Vde analysis (free parameters in bold). 


Symbol 

Definition 

Value 

Tf 

Mission duration [min] 

30 

N • 

sims 

Number of simulation runs 

10 

K 

Number of vehicles 

1 

Nt 

Number of time nodes 

50 

No, 

Number of parameter nodes 

25x25 

DU 

Canonical distance [m] 

100 

TU 

Canonical time [s] 

100 

UxV 

Vehicle dynamics file 

user-selected 

f 

Sonar operating frequency [kHz] 

[200, 450, 900} 

FOM 

Figure of Merit [dB] 

66 (median value) 

A 

Poisson Scan rate [scans/s] 

0.5 (median value) 

cr 

Signal excess Pd uncertainty [dB] 

9dB 

Heov 

Horizontal field of view [deg] 

90 (minimum value) 

Vfov 

Vertical field of view [deg] 

10 (median value) 

Vde 

Vertical mounting angle [deg] 

[- 25 : 5 :- 5 ] 


Additional simulations were condueted to determine the best mounting angles for different 
FLS deployed from a REMUS 100 AUV operating at 3 meters above the sea floor. An 
analysis of these results determined that the optimal mounting angles for detecting bottom 
mines from a survey altitude of 3 meters with the 450 kHz and 900 kHz FLS are -3 degrees 
and -5 degrees, respectively. The 200 kHz FLS was exempted from this analysis since the 
REMUS 100 AUV is unable to deploy a sensor of this size. The RID missions presented in 
Section 4.1 utilized these values for Vde- 

Figure 5.6 reveals that simulated survey missions with the 200 kHz FLS at down-tilt angles 
less than nine degrees produced only a handful of feasible trajectories. Although this may 
explain the relatively large standard deviation in with these mounting angles (see, e.g., 
the error bars for Vde = -6 and Vde = -8), more analysis is needed, perhaps using fewer 
time nodes, to verify this result. This is especially important since Vde has such a large 
impact on search performance. 
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MCM Risk vs. Sonar Vertical Mounting Angie (V^^) 
(1 USV searching 4 km^ area) 



Figure 5.5. Average single-vehicle search performance vs. mounting angle 

Vde- 


Fraction of Feasible Trajectories from 10 Vertical Mounting Angle (\^^) Simulations 



Figure 5.6. Fraction of single-vehicle Vde simulations with feasible trajecto¬ 
ries. 
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5.3.2 Figure of Merit 

For each of the three FLS considered in this analysis, Vde was set to the optimal value 
determined in Section 5.3.1. The other parameter values are listed in Table 5.3, where the 
bold array notation indicates the range of the independent variable under investigation, i.e., 
simulations are performed for FOM = 48,51,54,... ,75 dB. For each of the FOM values 
listed in Table 5.3, Figure 5.7 plots the mean and standard deviation (depicted by solid 
lines and error bars, respectively) for the objective values from ten simulations. Figure 5.8 
shows the fraction of these search trajectories which meet our feasibility criteria. Search 
performance for the solver output Jout and the ODE-propagated solution Jode (recalculated 
on our 500 x 25 x 25 discretization baseline) are both provided for comparison purposes. 

Table 5.3. Simulation parameters for FOM analysis (free parameters in bold). 


Symbol 

Definition 

Value 

Tf 

Mission duration [min] 

30 

N sims 

Number of simulation runs 

10 

N, 

Number of vehicles 

1 

Nt 

Number of time nodes 

50 

No. 

Number of parameter nodes 

25x25 

DU 

Canonical distance [m] 

100 

TU 

Canonical time [s] 

100 

UxV 

Vehicle dynamics file 

user-selected 

f 

Sonar operating frequency [kHz] 

[200, 450, 900} 

FOM 

Figure of Merit [dB] 

[ 48 : 3 : 75 ] 

A 

Poisson Scan rate [scans/s] 

0.5 (median value) 

cr 

Signal excess Pd uncertainty [dB] 

9dB 

Hpov 

Horizontal field of view [deg] 

90 (minimum value) 

Vpov 

Vertical field of view [deg] 

10 (median value) 

Vde 

Vertical mounting angle [deg] 

(-6,-7,-11} 
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MCM Risk vs. Sonar Figure of Merit (FOM) 
(1 USV searching 4 km^ area) 



Figure 5.7. Average single-vehicle search performance vs. figure of merit 
(FOM). 


Fraction of Feasible Trajectories from 10 FOM Simulations 


1 

S? 0.8 

s 0.6 

(ji 

ro 

0.4 

0.2 


^■200 kHz 
I ^450 kHz 
^■900 kHz 


60 65 

Figure of Merit (FOM) [dB] 


Figure 5.8. Fraction of single-vehicle FOM simulations with feasible trajec¬ 
tories. 


As expected, analysis of the FOM parameter for each sonar did not reveal a single optimal 
value; bigger is better. Since FOM represents the maximum (noise-limited) range at which 
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echoes are detectable, it makes sense that search performance improves linearly with in¬ 
creasing FOM, especially since propagation losses expressed in dB become linear at long 
ranges (see Figure 3.1). We observe, however, that the slope of these curves are steeper 
for the longer-range sonars. This is due to their small mounting angles, which allow the 
sonar beam to cover more ground before reaching the sea floor. Sonar beams emitted at 
larger mounting angles reach the sea floor before the sonar can utilize all of the additional 
detection range produced by the higher FOM . 

5.3.3 Poisson Scan Rate 

For each of the three FLS considered in this analysis, Vde was set to the optimal value 
determined in Section 5.3.1. The other parameter values are listed in Table 5.4, where the 
bold array notation indicates the range of the independent variable under investigation, i.e., 
simulations are performed for X = 0.1,0.2,0.3,..., 1.0 scans/second. For each of the X val¬ 
ues listed in Table 5.4, Figure 5.9 plots the mean and standard deviation (depicted by solid 
lines and error bars, respectively) for the objective values from ten simulations. Figure 5.10 
shows the fraction of these search trajectories which meet our feasibility criteria. Search 
performance for the solver output Jout and the ODE-propagated solution Jode (recalculated 
on our 500 x 25 x 25 discretization baseline) are both provided for comparison purposes. 


119 



Table 5.4. Simulation parameters for A analysis 


(free parameters in bold). 


Symbol 

Definition 

Value 

Tf 

Mission duration [min] 

30 

N • 

^^sims 

Number of simulation runs 

10 

K 

Number of vehicles 

1 

Nt 

Number of time nodes 

50 

No, 

Number of parameter nodes 

25x25 

DU 

Canonical distance [m] 

100 

TU 

Canonical time [s] 

100 

UxV 

Vehicle dynamics file 

user-selected 

f 

Sonar operating frequency [kHz] 

[200, 450, 900} 

FOM 

Figure of Merit [dB] 

66 (median value) 

A 

Poisson Scan rate [scans/s] 

[O.lrO.lrl.O] 

cr 

Signal excess Pd uncertainty [dB] 

9dB 

Hpov 

Horizontal field of view [deg] 

90 (minimum value) 

ypov 

Vertical field of view [deg] 

10 (median value) 

Vde 

Vertical mounting angle [deg] 

(-6,-7,-11} 


Figure 5.9 reveals the benefit of increasing the number of scans per unit time, but this im¬ 
pact is not as large as expected. Search performance increases rapidly (i.e., Pnd decreases) 
as the Poisson Scan rate A is increased from its lowest value of 0.1 scans/second to a value 
of 0.5 scans/second; beyond this value, the benefits of increasing A begin to diminish. Note 
that these benefits are greater for the 200 kHz sonar than the higher-frequency models. 
Figure 5.10 is concerning, however, as more than half of the trajectories produced by this 
sonar (using 50 time nodes) are infeasible. Additional simulations with fewer time nodes 
are needed to ascertain the relationship, if any, between higher A values and infeasible 
search trajectories for the 200 kHz sonar model. 
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MCM Risk vs. Sonar Poisson Scan Rate 
(1 USV searching 4 km^ area) 



Figure 5.9. Average single-vehicle search performance vs. Poission rate A. 


Fraction of Feasibie Trajectories from 10X, Simulations 



Figure 5.10. Fraction of single-vehicle A simulations with feasible trajecto¬ 
ries. 
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5.3.4 Horizontal Field of View 

For each of the three FLS considered in this analysis, Vde was set to the optimal value 
determined in Section 5.3.1. The other parameter values are listed in Table 5.5, where the 
bold array notation indicates the range of the independent variable under investigation, i.e., 
simulations are performed for Hfov = 30,45,60,..., 165 degrees. For each of the Hpov 
values listed. Figure 5.11 plots the mean and standard deviation (depicted by solid lines 
and error bars, respectively) for the objective values from ten simulations. Figure 5.12 
shows the fraction of these search trajectories which meet our feasibility criteria. Search 
performance for the solver output Jout and the ODE-propagated solution Jode (recalculated 
on our 500 x 25 x 25 discretization baseline) are both provided for comparison purposes. 

Table 5.5. Simulation parameters for Hpov analysis (free parameters in 

bold). 


Symbol 

Definition 

Value 

Tf 

Mission duration [min] 

30 

N sims 

Number of simulation runs 

10 

N, 

Number of vehicles 

1 

Nt 

Number of time nodes 

50 


Number of parameter nodes 

25x25 

DU 

Canonical distance [m] 

100 

TU 

Canonical time [s] 

100 

UxV 

Vehicle dynamics file 

user-selected 

f 

Sonar operating frequency [kHz] 

[200, 450, 900} 

FOM 

Figure of Merit [dB] 

66 (median value) 

A 

Poisson Scan rate [scans/s] 

0.5 (median value) 

cr 

Signal excess Pd uncertainty [dB] 

9dB 

Hpov 

Horizontal field of view [deg] 

[ 30 : 15 : 165 ] 

Vpov 

Vertical field of view [deg] 

10 (median value) 

Vde 

Vertical mounting angle [deg] 

(-6,-7,-11} 


Increasing Hpov has a dramatic impact on search performance, effectively increasing the 
sonar’s sweep width and producing faster area coverage rates during a constant velocity 
survey operation. These benefits can justify the cost of additional staves to increase the 
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horizontal FOV of a blazed array system. Not eoineidentally, maximizing horizontal FOV 
was a major driver in the design of the ATLAS sonar, which has Hpov >180 degrees. 


MCM Risk vs. Sonar Horizontal FOV 
{1 USV searching 4 km^ area) 



Figure 5.11. Average single-vehicle search performance vs. horizontal field 
of view [Hpov)- 


Fraction of Feasible Trajectories from 10 HFOV Simulations 



Figure 5.12. Fraction of single-vehicle Hpov simulations with feasible tra¬ 
jectories. 
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5.4 Multi-Vehicle Search Performance vs. Team Compo¬ 
sition 

As we have already seen in Seetion 4.2 it may not be possible to aehieve a desired risk 
threshold with a single seareh vehiele. If the assets are available, a foree eommander ean 
simply inerease the number of searehers operating in an area. While this solution is straight¬ 
forward, it may be sub-optimal. Worse yet, this approaeh ean pull needed assets out of 
another area, slowing the overall MCM operation. Pre-mission analysis ean help identify 
the vehiele and sensor eharaeteristies that produee the best system eonfiguration for a given 
MCM seenario. The GenOC framework ean serve as a pre-mission analysis tool to support 
these kinds of trade studies. Moreover, the ability to ineorporate realistie vehiele and sonar 
models to optimize mission-speeifie seareh objeetives ean produee more information than 
planning tools based solely on eoverage rates. In this seetion, we demonstrate our frame¬ 
work’s ability to analyze the mission effeetiveness of different autonomous vehiele teams 
eondueting the mine deteetion survey deseribed in Seetion 4.2. We begin with the simple 
ease of a team eomprised of one or more identieal seareh assets, eaeh equipped with the 
same forward-looking sonar. 

5.4.1 Number of Searchers 

This analysis eompares the seareh performanee for a team of identieal SeaFox US Vs, all of 
whieh are equipped with one of the three different FLS models listed in Table 2.1, where 
Vde has been set to the optimal value determined in Seetion 5.3.1. The other parameter 
values are listed in Table 5.6, where the bold array notation indieates the range of the inde¬ 
pendent variable under investigation, i.e., simulations are performed for = 1,2,3,... ,5 
seareh vehieles. For eaeh team eonfiguration listed in Table 5.6, Figure 5.13 depiets the 
average seareh performanee, while Figure 5.14 shows the number of feasible trajeetories 
produeed after ten simulations. 

Based on this analysis, me ean make the following observations: 

• A team of two US Vs equipped with the same 450 kHz blazed array FLS outperform 
a single USV equipped with a 200 kHz eylindrieal array FLS. This eould realize 
substantial eost savings, when a eommereial off-the-shelf sonar (e.g., a BlueView 
P450-90) is eompared against a developmental sonar like the ATLAS. 
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Table 5.6. Simulation parameters for multi-USV analysis (free parameters in 
bold). 


Symbol 

Definition 

Value 

Tf 

Mission duration [min] 

30 

N ■ 

sims 

Number of simulation runs 

10 

K 

Number of vehicles 

[ 1 : 1 : 5 ] 

Nt 

Number of time nodes 

50 

No. 

Number of parameter nodes 

25x25 

DU 

Canonical distance [m] 

100 

TU 

Canonical time [s] 

100 

UxV 

Vehicle dynamics file 

user-selected 

FLS 

Sonar parameters file 

user-selected 


• Two 200 kHz FLS-equipped searchers are required for complete coverage in the 
time available, but three 450 kHz FLS-equipped searchers achieve the same search 
performance. Again, a three-vehicle team using commercial P450-90 sonar may be 
less expensive to operate than a two-vehicle team using ATLAS-like sonar. 

• The 900 kHz sonar performs poorly when mounted on a USV operating in waters 
this deep. This high-resolution sonar should only be employed by AUVs conducting 
RID missions, e.g., in Section 4.1. 

• Figure 5.14 indicates that optimal solutions for all of the configurations tested pro¬ 
duce at least one infeasible trajectory, and the problem is exacerbated as the number 
of vehicles increase. 

• It is important to remember that each data point in a plot like this corresponds to a 
number of optimal trajectories which can be quickly exploited when needed. This 
feature of our proposed methodology is emphasized in Figure 5.1. 
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MCM Risk vs. Number of Vehicles 
{USV(s) searching 4 km^ area) 



Figure 5.13. Average search performance vs. the number of vehicles on a 
team. 


Fraction of Feasible Trajectories from 10 Multi-USV Simulations 



Number of Vehicles 


Figure 5.14. Fraction of multi-vehicle simulations that produce feasible tra¬ 
jectories. 


126 


































CHAPTER 6: 

Conclusions and Future Work 


6.1 Conclusions 

Autonomous vehicles will play an increasingly important role in military operations for 
the foreseeable future, particularly when danger to personnel is greatest. Mine counter¬ 
measures (MCM) is one critical mission area in which these systems have already proven 
themselves. Driven by the need to detect, identify, and neutralize the threat of underwater 
mines—in a variety of challenging environments—the U.S. Navy has invested heavily in 
vehicle and sensor technologies for different MCM missions. This has produced a spectrum 
of complementary capabilities, most of which have (necessarily) been built around dedi¬ 
cated vehicle platforms or sonar systems. One thrust of current naval research is focused 
on combining these capabilities via collaborative teams of autonomous vehicles. The goal 
of this approach is to maximize overall mission effectiveness while overcoming individual 
system limitations. 

Motion planning algorithms for these heterogeneous vehicle teams must consider the capa¬ 
bilities and limitations of each team member. Since many MCM missions involve searching 
for mines with a specific sonar payload, the current practice is to 1) partition the search area 
into separate zones according to sensor capabilities, and 2) assign specific vehicle/sonar 
assets to methodically search each zone using traditional geometric, e.g., lawnmower, cov¬ 
erage patterns. Planning multi-vehicle missions in this manner resembles a scheduling 
problem in which individual vehicle plans are sequenced and executed according to estab¬ 
lished area coverage rates. Given sufficient time to conduct a search, these methods do 
guarantee complete coverage. They are overly conservative, however, and do not exploit 
the full potential of highly-mobile vehicles to shorten operational timelines or maximize 
sensor performance. Motivated by new numerical techniques for solving optimal control 
problems with parameter uncertainty, this dissertation applies these techniques to address 
sensor-based motion planning for MCM search missions by autonomous vehicle teams. 

Optimal control provides a mathematical framework for solving motion planning problems 
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with dynamic constraints and different performance eriteria. Reeent advances in numeri¬ 
cal solution techniques have made it possible to address parameterized uneertainty in the 
objeetive function of an optimal control problem, permitting the direet optimization of 
seareh objeetives using GenOC. This approaeh can generate feasible vehicle- and sensor- 
specific motion trajectories that maximize mine deteetion probabilities under operational 
constraints, ffowever, realistic vehicle and sensor models are essential to producing opera¬ 
tionally relevant results. 

Chapter 2 of this dissertation develops physies-based probabilistie deteetion models for 
various mine hunting sonars eurrently employed for MCM missions. In these models, de¬ 
teetion performance is not only a function of the sonar’s design criteria but also the three- 
dimensional trajeetory exeeuted by the search vehicle. Although simplifying assumptions 
are made to faeilitate rapid caleulations during trajeetory optimization, the simulated de¬ 
teetion performance of these models agrees with expectations from actual sonar systems. 
Therefore, they represent a good eompromise between the traditional deteetion models used 
in seareh theory, and full-seale aeoustie simulators used for sonar design and performanee 
prediction. 

In Chapter 3, we deseribe how various MCM missions can be formulated as optimal seareh 
problems in the GenOC framework. Chapter 4 highlights the flexibility of our modular 
approaeh by generating motion plans for MCM wide area survey and target reaequisition 
missions via different combinations of searcher, sensor, and target distribution models. Op¬ 
timal seareh teehniques are partieularly useful for creating sensor-based trajeetories whieh 
outperform eonventional eoverage patterns under resource or time constraints. Examples 
are provided in Chapter 4 that illustrate how optimal solutions ean be used to establish new 
performanee benehmarks for a given scenario. This makes GenOC an attraetive tool for 
pre-mission analysis as well. 

A novel eontribution of this dissertation applies GenOC to investigate inverse problems re¬ 
lating MCM mission performance to specifie sonar designs or employment methods. Chap¬ 
ter 5 describes an approaeh for rapidly solving a range of optimal search problems in this 
framework to support Monte Carlo analysis of alternative eonfigurations. This approaeh 
is used to determine the best diseretization level for aceurate, yet feasible solutions in a 
given seenario. Based on this result, additional Monte Carlo experiments are eonducted to 
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analyze search performance relative to various sonar design parameters, highlighting the 
utility of this approach. 

6.2 Recommendations for Future Work 

This dissertation represents the first application of GenOC to operationally relevant motion 
planning problems in mine countermeasures. The derivation of physics-based models that 
incorporate design characteristics and three-dimensional detection performance of actual 
sonar systems has opened the door to a wide array of important analysis questions, and 
this dissertation only scratches the surface. We have demonstrated the utility of the GenOC 
framework for analyzing several questions of immediate interest to the MCM community, 
but a number of avenues remain open for future work. 

First, additional analysis is needed to investigate test cases from Chapter 5 that produced 
relatively low fractions of feasible solutions. This was more prevalent in simulations which 
utilized the 200 kHz FLS model. Repeating these experiments with fewer time nodes may 
produce more feasible trajectories. Whereas variation was added to these simulations via 
the initial condition, additional experiments could inject random variation into the “shape” 
of an initial guess trajectory. This approach may produce more globally optimal solutions 
by exploring more of the state space during a given Monte Carlo simulation. Another 
technique worth investigating is a two-step process whereby an optimal trajectory from 
one simulation is used as the initial guess for a second simulation. We regard all of these 
suggestions as low-hanging fruit. At a minimum, however, increasing the total number 
of simulations is recommended in order to produce more statistically-significant analytic 
results. 

As mentioned in Section 4.3, the computational framework used for this dissertation is 
designed around solving fixed-time optimal control problems. While this problem formu¬ 
lation can address time-limited missions directly, it requires additional effort to analyze 
explicit questions about time, e.g., how long would it take a team of vehicles to achieve a 
risk threshold of 5%? Currently, this requires a Monte Carlo analysis of simulations with 
an array oiTf values to provide an approximate solution. Future work will assess the abil¬ 
ity of other optimal control software packages like DIDO to solve minimum-time versions 
of optimal search problems for MCM. This software can also facilitate a more detailed 
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sensitivity analysis of inverse problems using the approach referenced in Section 5.1. 

This dissertation has presented results highlighting some of the trade-offs associated with 
selecting a discretization scheme to numerically solve optimal search problems defined in 
continuous time and space. This has proved especially challenging for the sidescan sonar 
models developed in Chapter 2, which are characterized by a very high Poisson Scan rate 
(T) and a very narrow horizontal field of view {Hpov)- Initial simulation results with these 
sensor models require an exceedingly large number of time nodes to produce accurate area 
coverage, and this can make numerical solutions intractable for complex motion planning 
problems. Since SSS is so widely used in underwater search applications, future work will 
continue to address this modeling challenge. 

An interesting avenue for future work is investigating ways to exploit the capabilities of 
AUVs and FLS to execute three-dimensional trajectories for mine detection. The solution 
framework readily supports vehicle dynamics with arbitrary degrees of freedom, and the 
sensor models developed in this dissertation already operate in three-dimensional space. 
Thus far, however, we have only considered search trajectories in a two-dimensional plane, 
in keeping with current MCM practices. For this reason, we consider the exploration of 
three-dimensional search trajectories to be high-risk, high reward, since a rigorous analysis 
of this approach could potentially improve MCM CONOPS—especially against moored or 
in-volume mines. 

The preceding recommendations can all be addressed with minimal changes by the current 
GenOC framework. The following suggestions would require significantly more effort, 
but would make important contributions. Perhaps most egregious to sonar practitioners, 
our current implementation does not address environmental impacts (beyond propagation 
losses) on sonar performance. Operational relevance can be greatly improved by incorpo¬ 
rating reverberation effects, sea floor characteristics, and sound speed profiles into the sonar 
detection models. A significant challenge, however, is the computational effort required by 
this additional complexity, since sonar detection calculations occur in the inner loop of the 
optimization algorithm. Again, we seek a balance between over-simplified “cookie cutter” 
detection models, and sonar performance simulations that take hours (or days!) to execute. 

Even with the simplified detection models used in this dissertation, current algorithms are 
not yet practical for real-time trajectory generation by vehicle autopilots. The current ap- 
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proach is better suited to pre-mission analysis in order to establish performance benchmarks 
and facilitate CONOPS development. Optimal motion plans generated by the solver are es¬ 
sentially open loop mission profiles that must be implemented separately by each vehicle. 
Future work to enable real-time computation and periodic re-planning of optimal search 
trajectories are key to transitioning this technology for use on actual MCM vehicles. 

Periodic re-planning would also allow vehicles to perform Bayesian updates on the un¬ 
derlying target distribution model, in response to detection events, and reallocate search 
effort more efficiently. This capability could be enhanced via wireless radio or acoustic 
communications between vehicles, allowing MCM teams to conduct survey missions and 
RID missions concurrently. This could realize significant time savings over current prac¬ 
tice, which conducts these MCM operations in sequential phases. It could also pave the 
way for mission plans which incorporate more complex vehicle interactions, e.g., the rapid 
transport of AUVs to and from an operating area by USV or helicopter to minimize the 
time wasted during slow-speed transits. 
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APPENDIX A: 
Objective Function Gradients 


A.l Objective Function Gradients 

The SNOPT numeric optimization package can solve an optimal search problem more 
quickly and reliably when it is supplied with analytic expressions for the objective func¬ 
tion’s gradients [99]. For our problem, we must derive gradients with respect to the 
state variables dJIdx and control inputs dJ/du for the objective function given by Equa¬ 
tion (3.11). We encode optimal search problems for SNOPT using a MATLAB interface 
developed by Claire Walton. This interface, Software for Parameter-uncertainty Optimal 
Control (SPOC), conveniently allows users to specify separate gradients for the inner, 
and outer, G(z) = e~^, components of a GenOC problem’s running 
cost in Equation (1.2) [79]. Since dGjdz has the trivial expression -e~^, and dyidu = 0 
for this problem, the following section derives dR/dx only. These gradients are based on 
components defined in Chapter 2 for the instantaneous detection rate of Equation (2.22). 
Note that all variables are assumed to be in canonical, non-dimensional form, but the over¬ 
line notation used in Section 3.5 is not used here. Eikewise, the searcher states’ dependence 
on time is omitted to simplify the notation. 

Recall that y(f,m) = Ap{x,co)Faix,co)Fi;{x,co)Fr{x). Each term in this equation except 


the constant Poisson Scan rate X depends on 

the state vector 

{x,y,il/,rY. Therefore, the 

gradient is 






dR dy 

dy 

dy 

dy 

dy 

(A.l) 

1 

1 

dx 

dy 


dr 
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We can apply the product rule to calculate each element of the objective function gradient: 


dy I dp dFa dF^ dFr 

/ = d -fF^F.Fr + p^F.Fr + pFa^Fr + pF^F,^ 
dx \dx dx dx dx 

dy I dp dFa dFg dFy 

^ = d -fFaF.Fr + P^F.Fr + pFa^Fr + pFaF,^ 
dy \dy dy dy dy 

dy I dp dFa dFf, dF^ 

dif/ \dif/ ^ dip ^ dip ^ ^ dip 

dy I dp dFa dF^ dF^ 

/ = d -fFaF.F, + p^F,F, + pFa^F, + pFaF,^ 
dr \ dr dr dr dr 


(A.2) 


The following sections derive the expressions for dpjdx, dFaldx, dF^/dx, and dFyjdx. 


A.2 Gradients for Instantaneous Probability of Detection 

(dpldx) 

Recall from Equation (2.3) that SE{x,6j) = FOM - PL (D(T,m)), where PL {D{x,6j)) 
is the propagation loss as a function of the distance between searcher and target, D{x,cd). 
Substituting into Equation (2.8) yields 


p(T,m) = O 


SE {x,(d) 


(T 


)=o( 


FOM -PL{D(x,c5)) 


cr 


O(^). 


(A.3) 


Eetting ^ = {FOM - PL)/a, we can rewrite the cumulative normal distribution function 
<I>(^) in terms of the erf function, which yields 




1 + erf 



We can now use the chain rule to derive dpjdx as 


dp dp df dPLdD 
dx df dPL dD dx 


(A.4) 


(A.5) 


Using the analytic expression for the derivative of the erf function, we compute dpjdf as 


dp 


_i£2 





FOM-PL 

cr 


(A.6) 
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Meanwhile, using our definition for we eompute: 


1 

dPL a 


(A.7) 


Using D = \\io - x\\ = y (cux - ■>c)^ + {cOy - y)^ + - z)^ = + dy^ + dz^ from 

Equation (3.22), and PL = 201ogio {DU) + 201ogio {D) + aD from Equation (3.26), we 
eompute dPL/dD and dD/dx as follows: 


dPL 

~dD' 


20 



) + 201ogio (D) + aD 


ln{D) 

20—^ + aD 
In(lO) 

20 


In(lO) 
1 


In(lO) \D 


ln(D) + aD 


+ a, and 


(A.8) 


^ ^ ^ {oJx - JC)2 + {ojy - y)2 + - z)2j 

^ 1_ 2(m^-.r)(-l) _ 

^ ^(cU;c - + {a^y - y)^ + {oiz - 

dx 
~ ~~D' 


(A.9) 


Einally, we substitute Equations (A.6) through (A.9) into Equation (A.5) to obtain 


dp 

dx 


dx 


crD V2^ 


20 


+ a 


dx 


cr 


20 + aDln(10) 


[ D2ln(10) 


1 L F-PL \^ 

■2V (T ) 


1 r F-PL\^ 
■2V o- ) 


In a similar manner, using dDjdy = -dy jD, we ean eompute dpjdy as 


dp dp d^ dPLdD 
dy d^ dPL dD dy 


dy 20 + aDln(10) 
cr^Pln . f)^ln(lO) 


1 r F-PL p 
e i\ <T ) _ 


Moreover, dp jdip = dpjdr = 0 sinee p{x,i5) is a funetion of distanee only. 


(A. 10) 


(A.ll) 
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A.3 Gradients for Azimuth Angle Shaping Function 

{dFaldx) 

Recall from Equation (2.15) that Fa{x,u) is based on the body-fixed azimuth angle ^a, 
which depends upon the states x, y, and ifj as shown in Equations (2.10) through (2.12), 
but is independent of turn rate r, i.e., dFa/dr = 0. We can rewrite Equation (2.15) as 


F - F + F - 1 


(A. 12) 


where: 


Far = -^-r = —, and 

1 - 1 - ePaictL - a) 

1 1 

f = _=_ 

1 + gPa(o! - au) Ajj ' 

Note that the body-fixed superscript b has been omitted from the azimuth angle a to sim¬ 
plify notation. We can now use the chain rule to derive dFa/dx as 


(A.13) 
(A. 14) 


dFa dFa dAi da dFa dAu da / dFa dAi dFa dAu\ da 
dx dAi da dx dAjj da dx da dAjj da ) dx" 


where: 


and: 


dFa _ 

1 

1 

(A.16) 

dAi 

[1-1- gPa 

(ttL - a) j 2 ’ 

dFa _ 

1 

1 

(A.17) 

dAu 

Al [1 + eP- 

,(a - au))'^’ 

II 

1 + gPAaL - _ 

-- 

(A.18) 

dAu d 

da da 

j^l + gPcrio: - Q;(7)j 


(A. 19) 


da 

dx 


d d d 

—atan2(*Jy, ^dx) = atan2("J^, "dx) = — arctan 
dx dx dx 



(A.20) 
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Letting tp = dyidx, we can compute the derivative of the inverse tangent function with 
respect to a: as 


— arctan(^) = 
ox 


1 ^ ] 

1 d(p j 

( dx^ \ 

1 

\1 +^2j 

' dx ' 

^dx^ + dy^ / 

' dx 


(A.21) 


Meanwhile, since dx = tOx - x, the derivative of cp with respect to x is 


dx dx \dxj dx^ dx^ 


(A.22) 


Substituting Equations (A.21) and (A.22) into Equation (A.20) yields 


da I dx^ \ dy dy 

dx \dx^ + dy^j dx^ dx^ + dy^ 


(A.23) 


Eetting EaL - and Ea^ = we substitute Equations (A. 16) through 

(A. 19), with Equation (A.23), into Equation (A. 15) and obtain 


dFc _ pady EgL _ Egjj 

dx ~ dx^ + dy^[^l+Eg^)^ {l+Eg^f 


(A.24) 


In a similar manner: 

dFg _ dFg dAi da ^ dFg dAg da _ / dFg dAi ^ dFg ^A[/ \ da lA 251 

dy dAi da dy dAu da dy ydAi da dAu da ) dy" 


^ = -;f-atan2(^Jy, ^dx) = -^atanlCdy, ’^dx) = arctan (^ | . (A.26) 

dy dy dy dy \dx j 


Once again, we let tp = dyjdx and compute the derivative of the inverse tangent function 
with respect to y as 


— arctan(^) = 
dy 


1 ^ ) 

1 dip j 

( dx^ \ 

1 dip 


' dy ' 

ydx^ + dy^ j 

' dy 


(A.27) 
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Meanwhile, sinee dy = ojy - y,the derivative of cp with respeet to y is 


^ ^ ^ ^ _J_ 

dy dy \ dx j dx dx 

Substituting Equations (A.27) and (A.28) into Equation (A.26) yields 


da j 

' dx^ \1 

f ^ 1 

1 dx 

dy ' 

^dx^ + dy^j ' 

dx f 

^ dx'^ -1- dy^ 


(A.28) 


(A.29) 


Using the same definitions for and Ea^, we substitute Equations (A. 16) through 
(A. 19), with Equation (A.29), into Equation (A.25) and obtain 


dFg, _ padx 

^CKL _ 

dy ~ dx^ + dy^[(^\+E^^f + E^^f 


(A.30) 


The derivation of dEaldij/ proeeeds similarly, as follows: 


dEa dEa dAi da dE^ dAjj da 
dif/ dAi da dij/ dAu da dtff \dAi^ da dAu da j difj 


dEa dAi dEa dAu \ da 


(A.31) 


However, note that d^ajdif/ 4^ d^ajdif/, and we must use the azimuth angle resolved in 
body-fixed eoordinates, whieh is given by 


d^a d u u d 

— = —atan2( dy, dx) = — aretan 


^dy 

^dx 


(A.32) 


Refereneing Eigure 2.8, we ean expose this funetion’s dependenee on if/ by letting 

^dy "dycos(ij/) - "dxsiniij/) 

^dx "Jxeos(^) -I- "dy sm(ij/) 

Using the ehain rule. Equation (A.32) beeomes 

d’^a d d dp 

^ aretan (/3(^)) = — aretan(/3) —. 
dip dip dp dip 


(A.33) 


(A.34) 
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The derivative of the inverse tangent function is 


Aarc,a„OT = ^ 
while the derivative of its argument jS is 


"dx^ + '^dy^ 


(A.35) 


d/3 d "dycos(i//) - "dxsm(il/) 
dif/ dtff "dxcos(i//) + "dysm(ij/) 

["dy cos(i//) - "dxsm(i//)] [- "dx sm(i//) + "dy cos(i//)] 
["dxcos{tff) + "dy sm{i//)]^ 

[- "dy - "J;ccos(i/r)] ["dxcos(i//) + ’^dy sm(i//)] 

["dxcos(d/) + "dy sm(ij/)]^ 

["dycos(i//) - "dx sm(i//)]^ + ["dx cos(if/) + "dy sm(i//)]^ 
["dxcos(i//) + "dy 

_ "dy^ cos^(tff) + "dx^ sin^(i//) + "dx^ cos^{i//) + "dy^ sm^(i//) 
["dxcos(d/) + "dy sm(ij/)]^ 

"dx^ [cos^(^) + sin^(i/^)] + "dy^ [cos^((/^) + sin^(^)j 

["dx cos(tff) + "dy sm(i//)]^ 

_ _ "dx^ + "dy^ 

["dxcos{i//) + "dy 


(A.36) 


When we combine Equation (A.35) and Equation (A.36) under the chain rule in Equa¬ 
tion (A.34), the partial derivative of azimuth angle with respect to heading angle if/ simpli¬ 
fies to 


d^a [^dxcosiifj) + sin(i/f)]' 
dij/ ’^dx^ + "dy^ 


"dx^ + "dy^ 

["dxcos(t//) + "dy sm(i//)]' 


= -1. (A.37) 


Since the z-axes of the inertial and body-fixed reference frames are parallel, azimuth angles 
in both coordinate systems can be related by = "a - ij/, as shown in Eigure 2.8. 
Erom this relation, d^a/dij/ = -1, which is consistent with Equation (A.37). The gradient 
dFa/dif/ in Equation (A.31) therefore simplifies to 


= -Pa 




(A.38) 


139 



A.4 Gradients for Elevation Angle Shaping Function 

{dFsIdx) 

Recall from equation Equation (2.16) that is based on the elevation angle s, which 

depends upon the states x and y as shown in Equation (2.14), but is independent of heading 
if/ and turn rate r, i.e., dF^/diJ/ = dF^ldr = 0. We can rewrite Equation (2.16) as 


Fs — Fe^ + Fejj 1 , 


(A.39) 


where: 


1 


1 


F = — 

1 + El 


, and 


1 


1 


F = — - 

l + ePXe-su) Eu' 


(A.40) 

(A.41) 


Note that elevation angle s is identical in both the inertial and body-fixed reference frames, 
so superscripts are omitted to simplify notation. We can now use the chain rule to derive 
dFs/dx as follows: 


OFe _ OFe dEi ds dpE dEu ds _ / OFe OEl dF^ dEu\ ds 
dx dEi ds dx dEu ds dx \dEL ds dEu ds j dx" 


(A.42) 


where: 


^ = -J-= 1 

dEu [1 + 

^ = -^ = 1 

dEu [1-I- 


dEu 
ds 
dEu _ 
ds 


ds 



1 ^ _p^^Ps{eL-e)^ 


ds d ( -dz \ 

- = -arctanj^-^==j. 


(A.43) 

(A.44) 

(A.45) 

(A.46) 

(A.47) 
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Letting ^ = -dz! ^jdx^ + dy^, we can compute the derivative of the inverse tangent func 
tion with respect to ;c as 


d / X / 1 \ / dx^ + dy^ \ (dx^ + dy^ 

— arctan(^) = [^ 772 ] ^ = \dx^ + dy^ + dz^] ^ " ( ~D^ 

Meanwhile, since dx = cOx - x, the derivative of ^ with respect to x is 

d^ d I -dz \ 1 Idxdz ^ dxdz 

\ ^jdx^ + dy^j ^ {dx'^ + dy'^)^^^ {dx^ + dy^)^^^ 

Substituting Equations (A.48) and (A.49) into Equation (AAV) yields 

ds Idx^ + dy^ 

dx \ 

Eetting and we substitute Equations (A.43) through 

(A.46), with Equation (A.50), into Equation (A.42) and obtain 


dxdz 


dxdz 


{dx^ + dy^)^^^ ^Jdx'^ + dy‘^ 


(A.50) 


dx 


(A.48) 


dF^ _ psdxdz _ ^eu 

(D2) ^dx^ + dy^ [(1 + ^.7" (1+^.7" 


(A.51) 


In a similar manner: 

dy dEi ds dy dEu ds dy xdEi ds dEjj ds ) dy' 


ds d I -dz \ 

— = —arctanj^-^==j 


(A.53) 


Once again,we let ^ = -dz! ^jdx'^ + dy^ and compute the derivative of the inverse tangent 
function with respect to y as 


d / \ I dx^ + dy^ \ d^ 

— arctan(0 = [77^] ^ = (^77777^) ^ 


I dx^ + dy^\ d^ 

( ) 77 - 


(A.54) 
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Meanwhile, sinee dy = ojy - y,the derivative of ^ with respeet to y is 


d ( -dz \ 1 2dydz ^ dydz 

^y\ ■\Jdx^ + dy^j {dx^ + dy^)^^^ {dx^ + dy'^)^^^ 

Substituting Equations (A.54) and (A.55) into Equation (A.53) yields 

de Idx^ + dy^ 

dy \ 

Using the same definitions for and we substitute Equations (A.43) through (A.46), 
with Equation (A.56), into Equation (A.52) and obtain 


dydz 


dydz 


{dx^ + dy^)^^^ 


(D^) 


(A.56) 


dFs _ psdydz ^EL _ ^Eg 

dy ~ (D2) ^dx^ + dy^ [{l + E,^f 


(AAV) 


A.5 Gradients for Turn Rate Shaping Function {dFr/dx) 

Reeall from Equation (2.20) that Er{x) depends only on the seareher vehiele’s turn 
rate r. Therefore, the gradients with respeet to all other state variables are zero, i.e., 
dEr/dx = dEr/dy = E^/dif/ = 0. The derivative of Er with respeet to r is simply 


dr dr 


(A.58) 
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A.6 Gradients for Instantaneous Detection Rate {dy/dx) 

We can now substitute the gradients derived in Sections A.2 through A.5 into the objective 
function gradients of Equation (A.2): 


dy 

dx 


dy 

dy 


dy 


dy 

dr 


( 




A 


dp dFa dFf; dF/r 

-^FaF^Fr + p^F^Fr + pF^—Fr + pF^F^^ 

dx dx dx /dx 


/ 


A 


A 


A 


dp dFa dFg dF 

-^FaFsFr + P^—FfrFr + pFa—Fr + pFaFg^ 

dy dy dy ^y 

^. 


‘J 


/ 




dFr 

dr 


, and 


Thus, the gradients of the inner component of the objective function’s running cost are 


dy 

dx 

dy 

dy 

dy 

dy 

dr 


'dp dFa dFe 

d -^FaFsFr + P^—FsFr + pFa — Fr , , 
, dx dx dx 


dp 


dFa 


dF, 


= A —FaFgFr + P^—FsFr + pFa—F, 
yy dy dy 

^ ’ and 

dF, 


= A pFaFs 


dr 


(A.59) 

(A.60) 

(A.61) 

(A.62) 
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APPENDIX B: 

Hamming Cluster Batch Scripts 


The Hamming supercomputer, a “hybrid cluster” of computing cores available via the NFS 
HPC Center [148] was used extensively during the course of this dissertation. Although 
our optimal control framework utilizes a sequential optimization algorithm that doesn’t 
lend itself to massive parallelization, several simulations can be launched via the network 
to run on multiple computing cores simultaneously. This capability was a key enabler for 
conducting numerous trade studies at once, greatly accelerating the analysis process. 

Hamming jobs are managed by Simple Linux Utility for Resource Management (SLURM), 
an open source “cluster management and job scheduling system” [150]. In order to use this 
resource effectively, a generalized MATLAB routine was developed to accept pertinent 
simulation parameters via SLURM environment variables. Next, batch scripts were cre¬ 
ated to execute this MATLAB routine on multiple computing nodes. These scripts leverage 
SLURM’s job array mechanism, whereby each job consists of multiple computing tasks, 
each indexed by a taskID value. The batch scripts specify a range of values for the param¬ 
eter of interest in a Monte Carlo simulation run, and convert these values into environment 
variables via the taskID. These environment variables are accessible by the MATLAB rou¬ 
tine, regardless of the computing node it has been assigned to run on. A separate input 
argument determines the number of simulations to run for a given parameter configuration. 
Upon completion, SLURM sends an email notification with a listing of all generated output 
files. These files are then downloaded from Hamming via the network and subsequently 
analyzed in MATLAB. 
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B.l SLURM Batch Files 


B.1.1 Launch file-based simulations 

# !I hin !bash 

# 

# Calling syntax is 

# 

# sbatch —array =<indVar range> 

# -o $OUTFILE_PREFIX.o%Ayca -J 

# <8-char job name> multiUxVsonarFile . sbatch 

# 

# Request hamming resources... 

# 

#SBATCH — ntasks=l 

#SBATCH — cpu-freq = high 

#SBATCH — constraint=intel 

#SBATCH —time =72:00:00 

#SBATCH — mail-user = spkragel@nps . edu 

#SBATCH — mail —type =END 

export OMP_NUM_TBtREADS= 1 

echo " Args ..." 

echo outfilePrefix : ${OUTFILE_PREFIX} 

echo Nv: $ {SPOC_ARG_NV} 

echo Tf: $ {SPOC_ARG_TE} 

echo Nt: $ {SPOC_ARG_NT} 

echo Nw: ${SPOC_ARG_NW} 

echo Method: ${SPOC_ARG_METHOD} 

echo DU: ${SPOC_ARG_DU} 

echo TU: $ {SPOC_ARG_TU} 

echo UxVparams: $ {SPOC_ARG_PARAMS} 

echo SonarEile: ${SPOC_ARG_SONAR} 
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echo numSims: $ {SPOC_ARG_N_SIMS} 

if [[ -z "${SPOC_ARG_NV}" ]]; then 

SPOC_ARG_NV=$ {SLURM_ARRAY_TASK_ID} 
echo "Nv from task:" ${SPOC_ARG_NV} 
elif [[ -z "${SPOC_ARG_NT}" ]]; then 
SPOC_ARG_NT=$ {SLURM_ARRAY_TASK_ID} 
echo "Nt from task:" ${SPOC_ARG_NT} 
elif [[ -z "${SPOC_ARG_TF}" ]]; then 
SPOC_ARG_TF= 

$( echo " scale =0; ${SLURM_ARRAY_TASK_ID}* 10" | be) 
echo "Tf from task:" ${SPOC_ARG_TF} 

else 

echo "All variables provided by user." 

fi 


if [[ -z "${SPOC_ARG_N_SIMS}" ]]; then 
SPOC_SPOC_ARG_TFARG_N_SIMS = 1 

fi 

echo "Running" $ {SPOC_ARG_N_SIMS} "MATLAB sims per task." 

echo "Submitted as:" ${SLURM_JOB_NAME} "with ID:" 

$ {SLURM_ARRAY_JOB_ID }"_"${ SLURM_ARRAY_TASK_ID} 


# Load MATLAB module 
source /etc/profile 
module load app/matlab 

# 

# cd to directory where job was submitted from 
cd $SLURM_SUBMIT_DIR 

# Run MATLAB 

if [[ -n "${SPOC_ARG_TF}" ]]; then 


147 



echo "Matlab script started on" 
date 

matlab -singleCompThread -nodesktop -nosplash -nojvm 
"ed . . / GoodMultiRobotConstVelSealedMultiFLS ; pwd ; 

-r runHeadless (${SPOC_ARG_NV} ,${SPOC_ARG_TF} , 

$ {SPOC_ARG_NT} , $ {SPOC_ARG_NW} , 

’ $ {SPOC_ARG_METHOD} ’ , $ {SPOC_ARG_DU} , 

$ {SPOC_ARG_TU} , ’ $ {SPOC_ARG_PARAMS } ’ , 

’ $ {SPOC_ARG_SONAR} ’ , $ {SPOC_ARG_N_SIMS } , 

’ $ {OUTFILE_PREFIX } ’ , $ {SEURM_ARRAY_JOB_ID} , 
$ {SEURM_ARRAY_TASK_ID }); quit" 
echo "Matlab script completed on" 
date 

else 

echo "ERROR! runHeadless has missing argument (s )! " 


fi 


148 



B.1.2 Launch parameter-based simulations 

# !I hin !bash 

# 

# Calling syntax is 

# sbatch —array =<indVar range> 

# -o $OUTFILE_PREFIX.o%Ayca -J 

# <8—char job name> multiUxVsonar . sbatch 

# 

# Request hamming resources... 

# 

#SBATCH —ntasks =1 

#SBATCH — cpu-freq = high 

#SBATCH — constraint=inte I 

#SBATCH —time =12:00:00 

#SBATCH — mail-user = spkragel@nps . edu 

#SBATCH — mail —type =END 

export OMP_NUM_TBtREADS= 1 

echo " Args ..." 

echo outfilePrefix : ${OUTFILE_PREFIX} 


echo 

Nv: 

${SPOC_ARG_NV} 

echo 

Tf: 

${SPOC_ARG_TF} 

echo 

Nt: 

${SPOC_ARG_NT} 

echo 

Nw: 

${SPOC_ARG_NW} 


echo Method: ${SPOC_ARG_METHOD} 
echo DU: ${SPOC_ARG_DU} 
echo TU: $ {SPOC_ARG_TU} 
echo UxVparams: $ {SPOC_ARG_PARAMS} 
echo f: $ {SPOC_ARG_EREQ} 
echo FOM: $ {SPOC_ARG_FOM} 
echo Eamda: $ {SPOC_ARG_EAMDA} 
echo Sigma: ${SPOC_ARG_SIGMA} 
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echo Hfov: $ {SPOC_ARG_HFOV} 
echo Vfov: $ {SPOC_ARG_VFOV} 
echo Vde: $ {SPOC_ARG_VDE} 
echo numSims: $ {SPOC_ARG_N_SIMS} 

if [[ -z "${SPOC_ARG_FOM}" ]]; then 

SPOC_ARG_FOM=$ {SFURM_ARRAY_TASK_ID} 
echo "FOM from task:" ${SPOC_ARG_FOM} 
elif [[ -z "${SPOC_ARG_FAMDA}" ]]; then 
SPOC_ARG_FAMDA= 

$( echo " scale =0; ${SFURM_ARRAY_TASK_ID}* 10" | be) 
echo "FAMDA from task:" ${SPOC_ARG_FAMDA} 
elif [[ -z "${SPOC_ARG_HFOV}" ]]; then 
SPOC_ARG_HFOV=$ {SFURM_ARRAY_TASK_ID} 
echo "Hfov from task:" ${SPOC_ARG_HFOV} 
elif [[ -z "${SPOC_ARG_VDE}" ]]; then 
SPOC_ARG_VDE=-$ {SFURM_ARRAY_TASK_ID} 
echo "Vde from task:" ${SPOC_ARG_VDE} 

else 

echo "All variables provided by user." 

fi 


if [[ -z "${SPOC_ARG_N_SIMS}" ]]; then 
SPOC_ARG_N_SIMS = 1 

fi 

echo "Running" $ {SPOC_ARG_N_SIMS} "MATFAB sims per task." 

echo "Submitted as:" ${SFURM_JOB_NAME} "with ID:" 

$ {SFURM_ARRAY_JOB_ID }"_"${ SFURM_ARRAY_TASK_ID} 


# Load MATLAB module 
source /etc/profile 
module load app/matlab 


150 



# cd to directory where job was submitted from 
cd $SLURM_SUBMIT_DIR 

# Run MATLAB 

if [[ -n "${SPOC_ARG_TF}" ]]; then 
echo "Matlab script started on" 
date 

matlab -singleCompThread -nodesktop -nosplash -nojvm 
"cd . . / GoodMultiRobotConstVelScaledMultiFLS ; pwd; 

-r runHeadless (${SPOC_ARG_NV} ,${SPOC_ARG_TF} , 

$ {SPOC_ARG_NT} , $ {SPOC_ARG_NW} , 

’ $ {SPOC_ARG_METHOD} ’ , $ {SPOC_ARG_DU} , 

$ {SPOC_ARG_TU} , ’ $ {SPOC_ARG_PARAMS } ’ , 

$ {SPOC_ARG_FREQ} , $ {SPOC_ARG_FOM} , 

$ {SPOC_ARG_LAMDA} , $ {SPOC_ARG_SIGMA} , 

$ {SPOC_ARG_HFOV} , $ {SPOC_ARG_VFOV} , 

$ {SPOC_ARG_VDE} , $ {SPOC_ARG_N_SIMS } , 

’ $ {OUTFILE_PREFIX } ’ , $ {SEURM_ARRAY_JOB_ID} , 
$ {SEURM_ARRAY_TASK_ID }); quit" 
echo "Matlab script completed on" 
date 

else 

echo "ERROR! runHeadless has missing argument (s )! " 
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B.2 SLURM Job Arrays 


B.2.1 Mission Duration 

# !I hin !bash 

# 

# Master script for passing variables into slurm 

# batch scriptas environment variables 

# 

# Calling syntax is 

# launch_PvsTf. slurm <OUTFILE_PREFIX> <SonarFile> <N_SIMS> 

# 

# 

export OUTFILE_PREFIX=$l 

export SPOC_ARG_NV=l # number of UxVs 

export SPOC_ARG_TE= # final time [s] 

export SPOC_ARG_NT=50 # time nodes 

export SPOC_ARG_NW=25 # spatial nodes 

export SPOC_ARG_METHOD="PEE" # P) seudospectral or E)uler 

export SPOC_ARG_DU=100 # dist canonical unit [m] 

export SPOC_ARG_TU=100 # time canonical unit [s] 

export SPOC_ARG_PARAMS=" SeaEox" # file of dynamic params 
export SPOC_ARG_SONAR=$2 # sonar params file 

export SPOC_ARG_N_SIMS=$3 # number of sims to run 

sbatch —array =90-366:12 —output=$OUTEIEE_PREEIX. o%A_%a 
—job-name=TE$2 multiUxVsonarEile . sbatch 
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B.2.2 Time Discretization 

# !I hin !bash 

# 

# Master script for passing va 

# batch script as environment 

# 

# Calling syntax is 

# launch_PvsNT. slurm <OUTFILE 

# 

# 

export OUTFILE_PREFIX=$l 

export SPOC_ARG_NV=l 
export SPOC_ARG_TE=1800 
export SPOC_ARG_NT= 
export SPOC_ARG_NW=25 
export SPOC_ARG_METHOD=" PEE" 
export SPOC_ARG_DU=100 
export SPOC_ARG_TU=100 
export SPOC_ARG_PARAMS= "REMUS" 
export SPOC_ARG_SONAR=$2 
export SPOC_ARG_N_SIMS=$3 

sbatch —array =20-60:5 —outpu 
—job-name=NT$2 multiUx 


'iables into slurm 
’ariables 

PREFIX> <SonarFile > <N_SIMS> 


# number of UxVs 

# final time [ s ] 

# time nodes 

# spatial nodes 

# P) seudospectral or E)uler 

# diSt canonical unit [m] 

# time canonical unit [ s ] 

# file of dynamic params 

# sonar params file 

# number of sims to run 

=$OUTEIEE_PREEIX. o%A_%a 
/sonarEile . sbatch 
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B.2.3 Vertical Mounting Angle 

# !I hin !bash 

# 

# Master script for passing variables into slurm 

# batch script as environment variables 

# 

# Calling syntax is 

# launch PvsVDE. slurm <OUTFILE PREEIX> <F KHZ> <N SIMS> 


# 

# 

export OUTFILE_PREFIX=$l 

export SPOC_ARG_NV=l 
export SPOC_ARG_TE=1800 
export SPOC_ARG_NT=50 
export SPOC_ARG_NW=25 
export SPOC_ARG_METHOD=" PEE" 
export SPOC_ARG_DU=100 
export SPOC_ARG_TU=100 
export SPOC_ARG_PARAMS= "REMUS" 
export SPOC_ARG_EREQ=$2 
export SPOC_ARG_EOM=66 
export SPOC_ARG_EAMDA=0.5 
export SPOC_ARG_SIGMA=9 
export SPOC_ARG_HEOV=90 
export SPOC_ARG_VEOV=10 
export SPOC_ARG_VDE= 
export SPOC_ARG_N_SIMS=$3 


# number of UxVs 

# final time [ s ] 

# time nodes 

# spatial nodes 

# P) seudospectral or E)uler 

# diSt canonical unit [m] 

# time canonical unit [ s ] 

# file of dynamic params 

# sonar frequency [kHz] 

# figure of merit [dB] 

# Poisson rate [1/s] 

# uncertainty [dB] 

# horizontal FOV ]deg] 

# vertical FOV ]deg] 

# sonar mount angle ]deg] 

# number of sims to run 


sbatch —array =5-25:5 —output=$OUTEIEE_PREEIX . o%A_%a 
—job -name=PvVDE$2 multiUxVsonar . sbatch 
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B.2.4 Figure of Merit 

# !I hin !bash 

# 

# Master script for passing variables into slurm 

# batch script as environment variables 

# 

# Calling syntax is 

# launch_PvsFOM. slurm <OUTFILE_PREFIX> <F_KHZ> <N_SIMS> 

# 

# 

export OUTFILE_PREFIX=$l 

export SPOC_ARG_NV=l # number of UxVs 

export SPOC_ARG_TE=1800 # final time [s] 

export SPOC_ARG_NT=50 # time nodes 

export SPOC_ARG_NW=25 # spatial nodes 

export SPOC_ARG_METHOD="PEE " # P) seudospectral or E)uler 

export SPOC_ARG_DU=100 # dist canonical unit [m] 

export SPOC_ARG_TU=100 # time canonical unit [s] 

export SPOC_ARG_PARAMS=" SeaEox " # file of dynamic params 

export SPOC_ARG_EREQ=$2 # sonar frequency [kHz] 

export SPOC_ARG_EOM # figure of merit [dB] 

export SPOC_ARG_EA]VIDA=0.5 # Poisson rate [1/s] 

export SPOC_ARG_SIGMA=9 # uncertainty [dB] 

export SPOC_ARG_HEOV=90 # horizontal FOV ]deg] 

export SPOC_ARG_VEOV=10 # vertical FOV [deg] 

export SPOC_ARG_VDE=-l 1 # sonar mount angle ]deg] 

export SPOC_ARG_N_SIMS=$3 # number of sims to run 

sbatch —array =48-75:3 —output=$OUTEIEE_PREEIX . o%A_%a 
—job -name=PvEOM$2 multiUxVsonar . sbatch 
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B.2.5 Poisson Scan Rate 

# !I hin !bash 

# 

# Master script for passing variables into slurm 

# batch script as environment variables 

# 

# Calling syntax is 

# launch_PvsLAM. slurm <OUTFILE_PREFIX> <F_KHZ> <N_SIMS> 

# 

# 

export OUTFILE_PREFIX=$l 


export SPOC_ARG_NV=l 
export SPOC_ARG_TE=1800 
export SPOC_ARG_NT=50 
export SPOC_ARG_NW=25 
export SPOC_ARG_METHOD=" PEE" 
export SPOC_ARG_DU=100 
export SPOC_ARG_TU=100 
export SPOC_ARG_PARAMS="SeaEox" 
export SPOC_ARG_EREQ=$2 
export SPOC_ARG_EOM=66 
export SPOC_ARG_EAMDA 
export SPOC_ARG_SIGMA=9 
export SPOC_ARG_HEOV=90 
export SPOC_ARG_VEOV=10 
export SPOC_ARG_VDE= -11 
export SPOC_ARG_N_SIMS=$3 


# number of UxVs 

# final time [ s ] 

# time nodes 

# spatial nodes 

# P) seudospectral or E)uler 

# diSt canonical unit [m] 

# time canonical unit [ s ] 

# file of dynamic params 

# sonar frequency [kHz] 

# figure of merit [dB] 

# Poisson rate [1/s] 

# uncertainty [dB] 

# horizontal FOV ]deg] 

# vertical FOV ]deg] 

# sonar mount angle ]deg] 

# number of sims to run 


sbatch —array =1-10 —output=$OUTEIEE_PREEIX . o%A_%a 
—job -name=PvEAM$2 multiUxVsonar . sbatch 
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B.2.6 Horizontal Field of View 

# !I hin !bash 

# 

# Master script for passing variables into slurm 

# batch script as environment variables 

# 

# Calling syntax is 

# launch_PvsHFOV. slurm <OUTFILE_PREFIX> <F_KHZ> <N_SIMS> 

# 

# 

export OUTFILE_PREFIX=$l 


export SPOC_ARG_NV=l 
export SPOC_ARG_TE=1800 
export SPOC_ARG_NT=50 
export SPOC_ARG_NW=25 
export SPOC_ARG_METHOD=" PEE" 
export SPOC_ARG_DU=100 
export SPOC_ARG_TU=100 
export SPOC_ARG_PARAMS="SeaEox" 
export SPOC_ARG_EREQ=$2 
export SPOC_ARG_EOM=66 
export SPOC_ARG_EAMDA=0.5 
export SPOC_ARG_SIGMA=9 
export SPOC_ARG_HEOV 
export SPOC_ARG_VEOV=10 
export SPOC_ARG_VDE=-6 
export SPOC_ARG_N_SIMS=$3 


# number of UxVs 

# final time [ s ] 

# time nodes 

# spatial nodes 

# P) seudospectral or E)uler 

# diSt canonical unit [m] 

# time canonical unit [ s ] 

# file of dynamic params 

# sonar frequency [kHz] 

# figure of merit [dB] 

# Poisson rate [1/s] 

# uncertainty [dB] 

# horizontal FOV ]deg] 

# vertical FOV ]deg] 

# sonar mount angle ]deg] 

# number of sims to run 


sbatch —array =30-165:15 —output=$OUTEIEE_PREEIX . o%A_%a 
—job -name=PvHEOV$2 multiUxV sonar, sbatch 
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B.2.7 Number of Searchers 

# !I hin !bash 

# 

# Master script for passing variables into slurm 

# batch script as environment variables 

# 

# Calling syntax is 

# launch_PvsNV. slurm <OUTFILE_PREFIX> <SonarFile> <N_SIMS> 

# 

# 

export OUTFILE_PREFIX=$l 

export SPOC_ARG_NV # number of UxVs 

export SPOC_ARG_TE=1800 # final time [s] 

export SPOC_ARG_NT=50 # time nodes 

export SPOC_ARG_NW=25 # spatial nodes 

export SPOC_ARG_METHOD= " PPP " # P) seudospectral or E)uler 

export SPOC_ARG_DU=100 # dist canonical unit [m] 

export SPOC_ARG_TU=100 # time canonical unit [s] 

export SPOC_ARG_PARAMS=" SeaFox " # file of dynamic params 
export SPOC_ARG_SONAR=$2 # sonar params file 

export SPOC_ARG_N_SIMS=$3 # number of sims to run 

sbatch —array =1-5 —output=$OUTFIEE_PREFIX . o%A_%a 
—job-name=NV$2 multiUxVsonarFile . sbatch 
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B.2.8 Type of Searchers 

# !I hin !bash 

# 

# Master script for passing variables into slurm 

# batch script as environment variables 

# 

# Calling syntax is 

# launch_PvsNV. slurm <NV> <UxvTeamStr> <SonarTeamStr> 

# <N_SIMS> <OutFilePrefix> 

# 

# 

export OUTFILE_PREFIX=$5 


export SPOC_ARG_NV=$l 
export SPOC_ARG_TE=1800 
export SPOC_ARG_NT= 
export SPOC_ARG_NW=25 
export SPOC_ARG_METHOD=" PEE 
export SPOC_ARG_DU=100 
export SPOC_ARG_TU=100 
export SPOC_ARG_PARAMS=$2 
export SPOC_ARG_SONAR=$3 
export SPOC_ARG_N_SIMS=$4 

sbatch 


# number of UxVs 

# final time [ s ] 

# time nodes 

# spatial nodes 

# P) seudospectral or E)uler 

# diSt canonical unit [m] 

# time canonical unit [s] 

# file of dynamic params 

# sonar params file 

# number of sims to run 


— array=10-45:5 — output =$OUTEIEE_PREEIX . o%A_%a 
—job-name=$2$3 multiUxVsonarEile . sbatch 
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B.3 MATLAB Run Script 

function Results = runHeadless ( varargin ) 

9(% New 11/18/16 for Heterogeneous Teams 
% 7 mandatory inputs: 

% Nv: number of UxVs 

% (ex: ’SRR ’ for 3-vehicle team with 1 SeaFox and 2 REMUS) 
% Tf: final mission time [seconds] 

% Nt: number of time nodes 
% Nw: number of parameter nodes 
% Methodstring : 3 chars for TSS, 

% either P for PS, E for Euler (ex: PEE) 

% DU: distance canonical unit [in meters] 

% TU: time canonical unit [ in seconds [ 

% 

% vehicle params 
% EITHER: 

% Veh: characters denoting team composition , 

% S for SeaFox, R for Remus 

% OR: (need all 3; these are applied to all Nv vehicles) 

% V: constant velocity [m/s [ 

% T: Nomoto T [ s [ 

% K: Nomoto K [1/s] 

% 

% sonar params 
% EITHER: 

% Son: characters denoting sonar model 

% (ex: ’249’ for 200kHz on UxV#l, 

% 450kHz on UxV#2, 900kHz on UxV#3) 

% OR: (need all 7; these are applied to all Nv vehicles) 

% f: design frequency [kHz] 

% FOM: figure of merit [dB[ 

% LAMDA: Poisson rate [l/s[ 

% SIGMA: cum. detection uncertainty [dB[ 
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% HFOV: horizontal field of view [degrees] 

% VFOV: vertical field of view [degrees] 

% VDE: vertical mounting angle [degrees] 

% (negative is down) 

% 

% Nsims: number of simulations 

% <outfile prefix >: prepends results file name 
% jobID: the slurm array job ID 
% taskID: the slurm array element ID 
% 

% 

% This launches multiple SPOC simulations . 

% SPOC, by CLAIRE WALTON (2013). Some Rights Reserved. 
9 ^==================================================== 


for s = l:Nsims 

disp([’*** Running Simulation num2str(s), , ... 

numlstr (numSims ) , ’ ***’]); 

CONSTANTS. cpuInfo = cpuinfo ; 

CONSTANTS .Start. Time = now ; 

CONSTANTS. Start . Date = d at e s tr (CONSTANTS . S tar t . Time ); 

Results = SPOC( gMultiCvScaledMultiFLS ,... 

Discretization , Methods) 

CONSTANTS. Stop . Time = now ; 

CONSTANTS. Stop . Date = datestr (CONSTANTS. Stop . Time ); 

save ( filename , ’Results’); 

disp ([’Results file is: ’ filename]); 

end; 
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